mmcdara

3713 Reputation

17 Badges

5 years, 357 days

MaplePrimes Activity


These are answers submitted by mmcdara

A very simple way:

restart:
with(Statistics):
HX := Histogram(Sample(Normal(1, 2), 1000), color="Blue"):

# type 
# op(HX)
# to see how HX is constructed

PLOT(
  POLYGONS(                  # this is the "main" part of the structure
    op(op(1, HX)),           # Histograms build
    LEGEND("data set 1")     #.... this is an information you could add
  ),
  op(2..-1, HX)              # and this is the last part ot the Histogram structure
)

 

This could be a trick
(even if I know that @acer has already told me he doesn't like much variable names constructed this way)

restart

A:= Vector(3, i -> v__||i)

A := Vector(3, {(1) = `#msub(mi("v"),mi("1"))`, (2) = `#msub(mi("v"),mi("2"))`, (3) = `#msub(mi("v"),mi("3"))`})

(1)

v__1 := 2:
A

Vector[column]([[v__1], [v__2], [v__3]])

(2)

v__1 := 'v__1':
A

Vector[column]([[v__1], [v__2], [v__3]])

(3)

 

Download trick.mw

I agree with @acer : as the numbers in your expressions are floats and not integer, rational, algebraic or transcendent number, the natural choice (IMO) of a solution method is fsolve and not solve.

But you could also solve the formal equality a = (1 + r__x)^b wrt r__x and evaluate the result with the floating point values.

restart
sol := solve(a = (1 + r__x)^b, r__x);
                            /ln(a)\    
                         exp|-----| - 1
                            \  b  /    
eval(sol, [a=1.15, b=1.28858]) 
                          0.114562536


Complement:
As said above "as the numbers in your expressions are floats and not integer, rational..."
That suggests to do this 

rel := 1.15 = (1 + r__x)^1.28858:
ratrel := convert(rel, rational)
                  23             (64429/50000)
                  -- = (1 + r__x)             
                  20                          

and then to apply solve (which is, IMO, now the natural choice) to ratrel.
Unfortunately the computational time is abnormally long.

Consider these two slight variants of your problem, 

rel := 23/20 = (1 + r__x)^1.2890625:
ratrel := convert(rel, rational);
solve(ratrel)
                    23             (165/128)
                    -- = (1 + r__x)         
                    20                      
                 1    (128/165)   (37/165)    
                 -- 23          20         - 1
                 20                           
rel := 23/20 = (1 + r__x)^1.2861328125:
ratrel := convert(rel, rational);
solve(ratrel)
                   23             (1317/1024)
                   -- = (1 + r__x)           
                   20                        
                                            1024    
                 /     1317                \        
           RootOf\20 _Z     - 23, index = 1/     - 1

In the second case, solve is probably stucked because it tries to compute a 1024 th root of the unity.
In your case, after the conversion into rational, the expression Maple gets must be something like this

RootOf(20*_Z^64429-23, index = 1)^50000-1

allvalues(%) will return this same expression, but fsolve(%) gives 0.114562536.



What is very strange, is that solve does exactly the operation sol := solve(a = (1 + r__x)^b, r__x); I did in the beginning of my answer, but fails to go further.
See below:

infolevel[solve]:=4:
rel := 1.15 = (1 + r__x)^1.28858:
ratrel := convert(rel, rational):
solve(ratrel);
Main: Entering solver with 1 equation in 1 variable
Dispatch: dispatching to Radicals handler
Dispatch: radical substitution will create polynomials of very high degree in (1+r__x)^(64429/50000).  Use of (1+r__x)^symbol may be more efficient
Transformer:   solving for linear equation in r__x
Recurse: recursively solving 1 equations in 1 variables
Dispatch: handling polynomials of the form a*x^n-b

This looks to me like a weakness of the solve procedure.

Some points fixed

  • pn is a recursive procedure which returns the correct result iif n=1.
    For n > 1 use thisproc for further call to pn within pn (look to Maple's help pages).
     
  • It seems to me that 
    tmp := RHS(x)*pn(i, 1, t)
    is not correct and should be
    tmp := RHS(t)*pn(i, 1, t)
    

Once corrected this way there is no longer any errors (corrections in magenta).
But you have to check yourself if what I did is right.

restart; with(LinearAlgebra)

h1 := proc (x) options operator, arrow; piecewise(0 < x and x < 1/2, 1, 1/2 < x and x < 1, -1, 0) end proc;

proc (x) options operator, arrow; piecewise(0 < x and x < 1/2, 1, 1/2 < x and x < 1, -1, 0) end proc

(1)

hi := proc (j, k, t) local a, b, c, m; m := 2^j; a := k/m; b := (k+1/2)/m; c := (k+1)/m; return piecewise(a <= t and t < b, 1, b <= t and t < c, -1) end proc:

J := 2:

for i to N do for j to N do H[i, j] := eval(hd[i], t = T[j]) end do end do:

 

RHS := proc (x) options operator, arrow; piecewise(-2 <= x and x < -1, 1, -1 <= x and x < 0, -1, 0 <= x and x < 1/2, 2, 1/2 <= x and x < 1, -2, 0) end proc;

proc (x) options operator, arrow; piecewise(-2 <= x and x < -1, 1, -1 <= x and x < 0, -1, 0 <= x and x < 1/2, 2, 1/2 <= x and x < 1, -2, 0) end proc

(2)

R := Vector(N); TMP := Matrix(N, N); A := Matrix(N, N)

R := Vector(4, {(1) = 2., (2) = 2., (3) = -2., (4) = -2.})

 

TMP := Matrix(4, 4, {(1, 1) = 1/4, (1, 2) = 3/4, (1, 3) = -3/4, (1, 4) = -1/4, (2, 1) = 1/4, (2, 2) = 3/4, (2, 3) = -3/4, (2, 4) = -1/4, (3, 1) = 1/4, (3, 2) = 1/4, (3, 3) = 0, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = -1/4, (4, 4) = -1/4})

 

A := Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0})

(3)

for i to N do R[i] := evalf(RHS(T[i])); tmp := RHS(t)*pn(i, 1, t); for j to N do TMP[i, j] := eval(tmp, t = T[j]) end do end do:

 

A := Transpose(LinearSolve(Transpose(H+TMP), R)):

proc (t) options operator, arrow; HFloat(1.0)*piecewise(t <= 0, 0, t <= 1/2, (1/2)*t^2, t <= 1, t-(1/2)*t^2-1/4, 1 < t, 1/4) end proc

(4)

``

Download wqw1_mmcdara.mw

Your worksheet contains a lot of arrors, some that I, or anyone else, can fix, and others no one but you can correct.
Plus, there are some lacks, no BCs for instance.
Other point, expecting an analytical solution is very likely a dead end.
Finally, please stop using this 2D input mode and stop using unnecessary parentheses (some are misplaced in your worksheet).

Anyway see a first draft, the rest is up to you

restart:

N := 15

15

(1)

# do not write f(t) and f[i] at the same time


f := t -> sum(p^i*phi[i]*t, i = 0 .. N)

proc (t) options operator, arrow; sum(p^i*f[i]*t, i = 0 .. N) end proc

(2)

# I guess Re is the Reynomds number and doesn't represents the Re(al) part
# of the expression between parenthesis ????

de1:=
  (1-p)*(
    (1+c1)*diff(f(t), t $ 4)
    -n0*(1+c1)*diff(f(t), t $ 2)  
    -m0*diff(f(t), t $ 2)
  )
  +p*(
    (1+c1)*diff(f(t), t $ 4)
    -c1*diff(g(t), t $ 2)
    -Rey*(
      f(t)*diff(f(t), t $ 3)
     -diff(f(t), t $ 1)*diff(f(t), t$2)
    )
    -n0 *(1+c1)*diff(f(t), t $ 2)  
    -m0*diff(f(t), t $ 2)
)

-p*c1*(diff(diff(g(t), t), t))

(3)

# If I rewrote correctly de1 then its trivial solution is
dsolve(de1);

# please verify this

g(t) = _C1*t+_C2

(4)

de2 :=
  (1-p)*(
    c2 *diff(g(t), t $ 2)
    +c1*(diff(f(t), t $ 2)-2*g(t))
  )
  +p*(
    c2 *diff(g(t), t $ 2)
    +c1*(diff(f(t), t $ 2)-2*g(t))
    -c3*Rey*(f(t)*diff(g(t), t $ 1)
    -diff(f(t), t $ 1)*g(t))
  );

(1-p)*(c2*(diff(diff(g(t), t), t))-2*c1*g(t))+p*(c2*(diff(diff(g(t), t), t))-2*c1*g(t)-c3*Rey*((f[0]*t+p*f[1]*t+p^2*f[2]*t+p^3*f[3]*t+p^4*f[4]*t+p^5*f[5]*t+p^6*f[6]*t+p^7*f[7]*t+p^8*f[8]*t+p^9*f[9]*t+p^10*f[10]*t+p^11*f[11]*t+p^12*f[12]*t+p^13*f[13]*t+p^14*f[14]*t+p^15*f[15]*t)*(diff(g(t), t))-(f[0]+p*f[1]+p^2*f[2]+p^3*f[3]+p^4*f[4]+p^5*f[5]+p^6*f[6]+p^7*f[7]+p^8*f[8]+p^9*f[9]+p^10*f[10]+p^11*f[11]+p^12*f[12]+p^13*f[13]+p^14*f[14]+p^15*f[15])*g(t)))

(5)

local D:
de3 :=
  (1-p)*(
    diff(a(t), `$`(t, 2))
    +D*(diff(b(t), `$`(t, 2)))
  )
  +p*(
    diff(a(t), t$2)
    +D*diff(b(t), t$2)
    +Ph*diff(f(t), t)*a(t)
    -ph*f(t)*diff(a(t), t)
  )

(1-p)*(diff(diff(a(t), t), t)+D*(diff(diff(b(t), t), t)))+p*(diff(diff(a(t), t), t)+D*(diff(diff(b(t), t), t))+Ph*(f[0]+p*f[1]+p^2*f[2]+p^3*f[3]+p^4*f[4]+p^5*f[5]+p^6*f[6]+p^7*f[7]+p^8*f[8]+p^9*f[9]+p^10*f[10]+p^11*f[11]+p^12*f[12]+p^13*f[13]+p^14*f[14]+p^15*f[15])*a(t)-ph*(f[0]*t+p*f[1]*t+p^2*f[2]*t+p^3*f[3]*t+p^4*f[4]*t+p^5*f[5]*t+p^6*f[6]*t+p^7*f[7]*t+p^8*f[8]*t+p^9*f[9]*t+p^10*f[10]*t+p^11*f[11]*t+p^12*f[12]*t+p^13*f[13]*t+p^14*f[14]*t+p^15*f[15]*t)*(diff(a(t), t)))

(6)

de4 :=
  (1-p)*(
    diff(a(t), t$2)
    +D*diff(b(t), t$2)
  )
  +p*(
    diff(b(t), t$2)
    +S*diff(a(t), t$2)
    +Pm*diff(f(t), t)*b(t)
    -ph*f(t)*diff(b(t), t)
  )

(1-p)*(diff(diff(a(t), t), t)+D*(diff(diff(b(t), t), t)))+p*(diff(diff(b(t), t), t)+S*(diff(diff(a(t), t), t))+Pm*(f[0]+p*f[1]+p^2*f[2]+p^3*f[3]+p^4*f[4]+p^5*f[5]+p^6*f[6]+p^7*f[7]+p^8*f[8]+p^9*f[9]+p^10*f[10]+p^11*f[11]+p^12*f[12]+p^13*f[13]+p^14*f[14]+p^15*f[15])*b(t)-ph*(f[0]*t+p*f[1]*t+p^2*f[2]*t+p^3*f[3]*t+p^4*f[4]*t+p^5*f[5]*t+p^6*f[6]*t+p^7*f[7]*t+p^8*f[8]*t+p^9*f[9]*t+p^10*f[10]*t+p^11*f[11]*t+p^12*f[12]*t+p^13*f[13]*t+p^14*f[14]*t+p^15*f[15]*t)*(diff(b(t), t)))

(7)

sys1 := eval({de1, de2, de3, de4}, p = 1):

# Whether my expression of de1 is correct or not you have 4 odes and 3 unknowns

select(has, indets(sys1, function), diff)

{diff(a(t), t), diff(b(t), t), diff(diff(a(t), t), t), diff(diff(b(t), t), t), diff(diff(g(t), t), t), diff(g(t), t)}

(8)

# Correct this and come to us later with a correct set of BCs

# once done,...
#
# you will have a lot of parameters to set to numerical values for no
# formal solution will be found (likely)


params := convert(indets(sys1, name), list)

[D, Ph, Pm, Rey, S, c1, c2, c3, ph, t, f[0], f[1], f[2], f[3], f[4], f[5], f[6], f[7], f[8], f[9], f[10], f[11], f[12], f[13], f[14], f[15]]

(9)

# Two possibilities:
#
#  1/ set all these parameters to numerical values and apply dsolve
#     for instance

# data := [D=1, Ph=2, .....]
# sol := dsolve(eval({sys1} union BCs, data), numeric)

#  2/ let these parameters "free" and use a parameterized dsolve

# sol := dsolve({sys1} union BCs, numeric, parameters=params)

# for an explanation of way 2 please look here dsolve[numeric,interactive]

Download Hpm_try_19_mmcdara.mw


EDITED : your boundary conditions can,ot be written this way

cond[1][0] := f[0](-1) = -1, (D(f[0]))(-1) = 0, f[0](1) = 1, (D(f[0]))(1) = 0; g[0](-1) = 0, g[0](1) = 0, a[0](-1) = 1, a[0](1) = 1, b[0](-1) = 1, b[0](1) = 0

for neither  f[0], a[0], b[0], g[0] are functions of t.
Concerning a, b or g the correct writing must be of the form 

a(-1)=..., a(1)=..., ...g(1)=...

Concerning f, unless I did some mistakes in copying-pasting your odes into a 1D input mode worksheet, f(t) doesn(t appear in any of them. You then have 4 odes and only 3 unknowns (a, b and g).

restart:
with(combinat, cartprod):
TA := [A, B, C]:
TB := [X, Y, Z]:

cc := cartprod([TA, TB]):
Tournaments := NULL:
while not cc[finished] do 
  Tournaments := Tournaments, cc[nextvalue]() 
end do:
Tournaments
   [A, X], [A, Y], [A, Z], [B X], [B, Y], [B, Z],[C, X], [C, Y], [C, Z]

 

Equivalent Resistance only,

(I do not understand what 'Solve minimum diagonal surface distance of nodes (0,0,0)-(3,4,5). (sqrt(74))" means)

Here I use the package GraphTheory
This is a luxury approach but it's a generic one which can be applied on any other any electrical circuit "topography".
The main trick is to use a more general expression of Kirchoff's laws to build the desired eqyuations in an automatic way:

  • instead of writting, as usual, that "the difference of the electical potentials between two points A and B is a constant whatever the path taken from A to B"
  • I write that the total variation of the electrical potential along any cycle from the cycle basis of the graph is identically null.

Note your problem reduces to a full rank linear system of 12 equations.

Equivalent_resistance.mw

As "Maple experiences difficulties in loading the content of the file", here is the result

                                       156
                     R["equivalent"] = ---
                                       47 

I answer here your reply to Preben:

SENSITIVITY ANALYSIS

There exists a lot of way to do this depending on:
   1/ you want operate within a statistical or deterministic framework
   2/ you want to do what is named  a "local sensitivity analysys" (LSA) or  a "global sensitivity analysys" (GSA)
       (please look to wikipedia to have a quick understanding on what these terms mean)

The simplest LSA example every student knows resumes to compute first order partial differences at some point (E, T)
Each partial derivatives with respect to some parameter q gives an information of the sentitivity of the "output" Y (here Y = S_TEV - S_P)  wrt q.
Here diff(Y, q) is a measure of the sensitivity of Y wrt q at some chosen point (E=e, T=t). Of course it depends, a priori, on the point you choosed.

This same approach within a statistical framework consists in considering the parameters E and T as random variables with some distributions.
A simple case is "first order LSA": is m(E) and m(T) are the expectations of E and T,  V(E) and V(T) are theis variances Cov(T, E) their covatiances, 
then the first order variation of Y is
   V(Y) = diff(Y, E)^2*V(E) +  diff(Y, T)^2*V(T) + 2*Cov(E, T)*diff(Y, E)*diff(Y, T)
and, assuming that Y := (E, T) --> f(E, T),  the 1st order approcimation of the expectation of Y is
   m(Y) : f(m(E), m(T))
Then the Local Sensitivity Indices LSI(Y | q) of Y wrt q, assessed at point (E=m(E), T=m(T)) are:
   LSI(Y | E) = diff(Y, E)^2*V(E) / V(Y)
   LSI(Y | E) = diff(Y, T)^2*V(T) / V(Y)

This is basically what the Maple's package ScientificErrorAnalysis does.

For GSA can be done in many many different ways.
In the statistical framework the most known approach is to use a Monte-Carlo strategy (several others do exist).
In the deterministic framework GSA is based on a prior definition (not only one does exist)  of the "Total Variation" of a function. One of the most common used approaches is based on the computation of Sobol's sensitivity analysis.

This means that no reasonable person can give you advice on how to do a sensitivity analysis (SA): it is for you to define more precisely what you mean by SA, only then can I, or others, give you useful advice.

SA.mw

You use some version of Windows and you have probably copied-pasted the path wich appears in some window.
In this case the path is in "Windows format", that is with backslashes, and not in Unix/LinuxMax OSX format (withs slashes instead)).
The first thing to do is to change all backslashes into slashes: (C://Users/....)
 



The position of the legend and the line spacings can be adjusted as you want.
Hint : 

  • op~(2, [plottools:-getdata(p1)]) gets the x and y ranges of all the plots
  • op~(2, op~(2, [plottools:-getdata(p1)])) gets the y ranges of all the plots
     

restart:

with(plots):

equ1 := BesselJ(sqrt(17)/2, 10*sqrt(t)*sqrt(2))/t^(1/4) + BesselY(sqrt(17)/2, 10*sqrt(t)*sqrt(2))/t^(1/4):

equ2 := BesselJ(sqrt(17)/2, 10*sqrt(t)*sqrt(2))/t^(1/4) + 5*BesselY(sqrt(17)/2, 10*sqrt(t)*sqrt(2))/t^(1/4):

equ3 := BesselJ(sqrt(17)/2, 10*sqrt(t))/t^(1/4) + 5*BesselY(sqrt(17)/2, 10*sqrt(t))/t^(1/4):

tmax   := 30:
colors := ["Red", "Violet", "Blue"]:

p1 := plot([equ1, equ2, equ3], t = 0 .. tmax, labels = [t, T[2](t)], tickmarks = [0, 0], labelfont = [TIMES, ITALIC, 12], axes = boxed, color = colors):

ymin := min(op~(1, op~(2, op~(2, [plottools:-getdata(p1)])))):
ymax := max(op~(2, op~(2, op~(2, [plottools:-getdata(p1)])))):
dy   := 2*ymax:

legend1 := typeset(C[3] = 1, ` , `, C[4] = 1, ` , `, Omega^2 = 50):
legend2 := typeset(C[3] = 1, ` , `, C[4] = 5, ` , `, Omega^2 = 50):
legend3 := typeset(C[3] = 1, ` , `, C[4] = 5, ` , `, Omega^2 = 25):

p2 := seq(textplot([tmax-2, ymax-k*dy/20, legend||k], align=left), k=1..3):

p3 := seq(plot([[tmax-2, ymax-k*dy/20], [tmax-1, ymax-k*dy/20]], color=colors[k]), k=1..3):
display(p1, p2, p3, view=[default, -ymax..ymax], size=[800, 500])

 

 


 

Download Legend_Inside.mw



Leap year not accounted for (every february here has 28 days)
 

restart:
month := [seq(cat("0", k), k=1..9), seq(cat("1", k), k=0..2)]:
day   := convert~([$1..9, seq(seq(cat(i, k), k=1..9), i=1..2), 30, 31], string):
bounds := [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]:
alias(CR=combinat:-randperm):
m := CR(month)[1];
d := CR(day[1..parse(m)])[1];
y := parse(StringTools:-Reverse(cat(d, m)));
                              "09"
                              "7"
                              907
cat(cat(d, "/", m), "/", y)
                           "7/09/907"

 


For non linear fit I usually prefer using Optimization:-NLPSolve instead of Statistics:-NonlinearFit.
As a rule Optimization:-NLPSolve is much faster than Statistics:-NonlinearFit and, as it is very flexible to define constraints, it avoids doing this loop on n.

Using  Optimization:-NLPSolve will enable you to select easily the "correct best solutions" (in the second attempt I used c >=-100).

I let you compare the execution time using  Optimization:-NLPSolve or your code (in black at the end of the file).

 

restart;

with(Statistics):

j := Matrix(20, 3, [[70, 8, 8.468140006], [70, 10, 4.105515432], [70, 12, 2.36261199], [70, 14, 1.422093923], [70, 16, 0.9], [100, 8, 20.47249229], [100, 10, 9.618450629], [100, 12, 5.360869165], [100, 14, 3.399312905], [100, 16, 2.640399788], [130, 8, 35.90466304], [130, 10, 17.62958097], [130, 12, 9.828362586], [130, 14, 5.866863694], [130, 16, 3.799262645], [160, 8, 57.31814648], [160, 10, 34.49692774], [160, 12, 15.39340528], [160, 14, 9.991012951], [160, 16, 6.049343013]]):

# last command killed after 30 seconds
# eq := NonlinearFit(u*h^n*exp(-b*d)/(-c*d + h), j, [h, d], output = [parametervalues, residualsumofsquares]);

# define the model f(h, d) and the objective function obj (sum of squared residuals)

f   := (h, d) -> u*h^n*exp(-b*d)/(-c*d + h):
obj := add( (f(j[k, 1], j[k, 2]) - j[k, 3])^2, k=1..numelems(j[..,1]) ):

# Use Optimization:-NLPSolve instead of Statistics:-NonlnearFit (the former is more versatile)

opt := Optimization:-NLPSolve(obj, {n>=2, n<=3}, method=sqp, iterationlimit=1000)

[21.6509289224931862, [b = HFloat(0.2028991836842215), c = HFloat(-6477.6480267378065), n = HFloat(2.3544626835278906), u = HFloat(99.32223640618867)]]

(1)

BestModel := eval(f(h, d), opt[2])

HFloat(99.32223640618867)*h^HFloat(2.3544626835278906)*exp(-HFloat(0.2028991836842215)*d)/(HFloat(6477.6480267378065)*d+h)

(2)

plots:-display(
  ScatterPlot3D(j, symbol=solidsphere, color=green, symbolsize=15),
  plot3d(BestModel, h=(min..max)(j[..,1]), d=(min..max)(j[..,2]), style=wireframe, color=blue)
);

 

# add an extra constraint on c

opt := Optimization:-NLPSolve(obj, {n>=2, n<=3, c >=-100}, method=sqp, iterationlimit=1000);
BestModel := eval(f(h, d), opt[2]):
plots:-display(
  ScatterPlot3D(j, symbol=solidsphere, color=green, symbolsize=15),
  plot3d(BestModel, h=(min..max)(j[..,1]), d=(min..max)(j[..,2]), style=wireframe, color=blue)
);

[20.2914105165715100, [b = HFloat(0.24364537500713906), c = HFloat(-20.866765615558734), n = HFloat(2.7593674638195393), u = HFloat(0.11018019031068993)]]

 

 

restart;

with(Statistics):

j := Matrix(20, 3, [[70, 8, 8.468140006], [70, 10, 4.105515432], [70, 12, 2.36261199], [70, 14, 1.422093923], [70, 16, 0.9], [100, 8, 20.47249229], [100, 10, 9.618450629], [100, 12, 5.360869165], [100, 14, 3.399312905], [100, 16, 2.640399788], [130, 8, 35.90466304], [130, 10, 17.62958097], [130, 12, 9.828362586], [130, 14, 5.866863694], [130, 16, 3.799262645], [160, 8, 57.31814648], [160, 10, 34.49692774], [160, 12, 15.39340528], [160, 14, 9.991012951], [160, 16, 6.049343013]]):

minRSS := +infinity:
for n from 2 by 0.001 to 3 do

    eq := NonlinearFit(u*h^n*exp(-b*d)/(-c*d + h), j, [h, d], initialvalues = [b = 0.1, c = -30, u = 0.1], output = [parametervalues, residualsumofsquares]);

    if eq[2] < minRSS then

        #print(n, eq[2], eq[1]);
        optRSS := [''n''=n, ''minRSS'' = eq[2], eq[1]]:
        minRSS := eq[2]:

    end if;

end do:

optRSS;

[n = 2.759, minRSS = 20.29141110, [b = HFloat(0.24360946347455648), c = HFloat(-20.898174204286263), u = HFloat(0.11043917832727822)]]

(3)

 


 

Download ConstrainedNonLinearFit.mw

 

To complete @tomleslie 's answer

If your matrix is that big, save it in an m file (binary file): you will spare time and bytes.

Let M the name of the matrix and f := "//../MyFile.m" some m file

  • save M, f   will save the matrix

To read it in a new worksheet

  • define f := "//../MyFile.m"
  • read f  will read the content of f
    This command automatically creates a matrix whose name is the same than the name used in the save command.
     

To verify this :

restart:
f := "//..../MyFyle.m";
read f;          # no echo of what is done
anames(user);    # will display the user variables, here f and M
                 # Useful if you don't remember the name the matrix had in the initial worksheet
                 # or if you share your work with other people

 


Not perfect but I have only Maple 2015 right now.
(advice: browse the questions to seek the works @acer has already done on related problems)

Download Colors.mw


1 2 3 4 5 6 7 Last Page 1 of 30