tomleslie

13841 Reputation

20 Badges

15 years, 1 days

MaplePrimes Activity


These are replies submitted by tomleslie

that if I wanted to

Please, how can I get various pp for different xi, say xi = -2..5, and then plot pp versus xi?

I wouldn't start from here! The way your code is structured makes it difficult to achieve this. A crude but functional approach is just to do a "cut-and-paste" job embedding all of your code within a procedure, pass in a value of 'xi' as an argument and have it return the value you want.

Anythig other than this crude "cut-and-paste" would take longer than I am prepared to spend on this!

The attached appears to do what you want, and allows you to constrcuct a graph of

min( seq( Re(gamma1*alpha/(4*k__1*qq[i]^2/d^2-alpha3*xi/eta__1)), i=1..m) )

as 'xi' varies.

NB the final graph renders a lot better in the worksheet than it does on this site

  restart;
  doCalc:= proc( xi )
                 uses ArrayTools, Student:-Calculus1, LinearAlgebra,
                      ListTools, RootFinding, ListTools:
                 local gamma1:= .1093,
                       alpha3:= -0.1104e-2,
                       k__1:= 6*10^(-12),
                       d:= 0.2e-3,
                       theta0:= 0.1e-3,
                       eta__1:= 0.240e-1,
                       alpha:= 1-alpha3^2/(gamma1*eta__1),
                       c:= alpha3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha3*xi/eta__1)),
                       theta_init:= theta0*sin(Pi*z/d),
                       n:= 10,
                       g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,
                       qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,
                       soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,
                       AllAsymptotes;

                 g:= q-(1-alpha)*tan(q)-c*tan(q):
                 f:= subs(q = x+I*y, g):
                 b1:= evalc(Re(f)) = 0:
                 b2:= evalc(Im(f)) = 0:
                 qstar:= (fsolve(1/c = 0, q = 0 .. infinity)):
                 OddAsymptotes:= Vector[row]([seq(evalf(1/2*(2*j + 1)*Pi), j = 0 .. n)]);
                 ModifiedOddAsym:= abs(`-`~(OddAsymptotes, qstar));
                 qstarTemporary:= min(ModifiedOddAsym);
                 indexOfqstar2:= SearchAll(qstarTemporary, ModifiedOddAsym);
                 qstar2:= OddAsymptotes(indexOfqstar2);
                 AreThereComplexRoots:= type(true, 'truefalse');
                 try
                      soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = 0 .. infinity});
                      soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = -infinity .. 0});
                      qcomplex1:= subs(soln1, x+I*y);
                      qcomplex2:= subs(soln2, x+I*y);
                 catch:
                       AreThereComplexRoots:= type(FAIL, 'truefalse');
                 end try;
                 OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);
                 AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));
                 if   AreThereComplexRoots
                 then gg:= [ qcomplex1,
                             qcomplex2,
                             op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),
                             seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)
                          ];
                 elif not AreThereComplexRoots
                 then gg:= [ op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),
                             seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)
                           ];
                 end if:
                 qq:= MakeUnique(gg):
                 m:= numelems(qq):
                 return min( seq( Re(gamma1*alpha/(4*k__1*qq[i]^2/d^2-alpha3*xi/eta__1)), i=1..m) );

            end proc:

#
# Run the calculation for a couple of values of 'xi'
#
  doCalc(-0.01);
  doCalc(-2.5);

.2045951722

 

-2.798195029

(1)

#
# Plot results of calculation for a range of xi-values.
# This takes a minute or so on my machine
#
  plots:-pointplot([seq( [j, doCalc(j)], j=-2.5..0, 0.1)], style=pointline);

 

 

Download cutPaste.mw

 

just use

pp:= min( seq( Re(gamma1*alpha/(4*k[1]*qq[i]^2/d^2-alpha3*xi/eta[1])), i=1..m) );

??

@pik1432 

I now understand what you are trying to achieve. You want to construct a differential equation in terms of Linear Differential Operators. Maple has a whole set of commands for doing this. See the help page at ?DEtools, in particular the section headed "Commands for Differential Operators".

You may also want to read the help page at ?Linear Differential Operators.

I'm not convinced how much this will help with your original problem since it appears to have two dependent variables, Delta__delta(t) and Delta__psi__d(t)

The attached shows some simple uses and manipulations as examples

  restart;
  with(DEtools):
  _Envdiffopdomain:=[Dt,t]:
#
# A good idea would be to start by reading the help page
# on Linear differential operators, here
#
? Linear Differential Operators;

#
# Define a couple of such operators, and construct a third
# by multiplication
#
  L1:=a__0+a__1*Dt;
  L2:=b__0*t+b__2*t^3*Dt^2;
  L3:=mult(L1, L2);

Dt*a__1+a__0

 

Dt^2*b__2*t^3+b__0*t

 

a__1*b__2*t^3*Dt^3+a__0*b__2*t^3*Dt^2+3*a__1*b__2*t^2*Dt^2+a__1*b__0*t*Dt+a__0*b__0*t+a__1*b__0

(1)

#
# Tidy the above expression by using collect() on Dt
#
  L4:=collect(L3, Dt);

a__1*b__2*t^3*Dt^3+(a__0*b__2*t^3+3*a__1*b__2*t^2)*Dt^2+a__1*b__0*t*Dt+a__0*b__0*t+a__1*b__0

(2)

#
# Produce the factors of L4. Note that each factor is
# also a linear differential operator
#
  allvalues~(DFactor(L4));

[a__1*b__2*t^3*Dt+b__2*t^2*(a__0*t+3*a__1), Dt+(1/2)*(2*(-(1/4)*(4*b__0-b__2)/b__2)^(1/2)+1)/t, Dt+(1/2)*(-2*(-(1/4)*(4*b__0-b__2)/b__2)^(1/2)+1)/t, Dt-(1/2)*(2*(-(1/4)*(4*b__0-b__2)/b__2)^(1/2)+1)/t, Dt-(1/2)*(-2*(-(1/4)*(4*b__0-b__2)/b__2)^(1/2)+1)/t]

(3)

#
# Apply the differential operator L4 to the dependent variable
# say f(t), and convert to a differential equation
#
  de1:=diffop2de(L4, f(t));

(a__0*b__0*t+a__1*b__0)*f(t)+a__1*b__0*t*(diff(f(t), t))+(a__0*b__2*t^3+3*a__1*b__2*t^2)*(diff(diff(f(t), t), t))+a__1*b__2*t^3*(diff(diff(diff(f(t), t), t), t))

(4)

#
# Alternatively generate the differential operator from
# the above differential equation - whihc should get back
# to something like L4 above
#
  L5:= de2diffop( de1, f(t) );

a__1*b__2*t^3*Dt^3+(a__0*b__2*t^3+3*a__1*b__2*t^2)*Dt^2+a__1*b__0*t*Dt+a__0*b__0*t+a__1*b__0

(5)

 

Download LDOs.mw

@pik1432 

You state

This question is a step further from the previous one. What I was trying to ask here is how to get the expression, presented by 'desired' in the worksheet attached to my question. 

You cannot leave an operator like 'D' dangling with nothing to operate on - it doesn't make sense. If you really need a trivial illustration then consider

x+sin(x)

whihch apparently you wish to be able to write as

x*(1+sin());

Do you really expect this to work? I mean really!

@alfarunner 

can have only one index. But there is nothing to stop you "nesting" sequences, whihc is what I didd in my original answer (and what Rouben did in his answer)

@Rouben Rostamian  

I'd call that coming up with your own three spatial boundary conditions - which may (or may not!) be the three spatial boundary conditions whihc the OP requires!

@Carl Love 

kernelopts(asserlevel) returns '0' - the default, so I'm pretty sure my problem is not associated with "extra" assertion checking.

I have come to the conclusion that it actually is an OS issue. I'm still(!!) running Windows 7. Technically,the system requirements for both Maple 2021 and Maple 2022, specify Windows 10 (although, so far I've never had a problem running these on Windows 7).

The DEQueue() object was only introduced in Maple 2021 - so maybe the associated code was never tested on a Win7 machine?

BTW On my machine, some commands associated  with DEQueue() "work" and some don't. Other than the fact that problems only *seem* to occur when dealing with zero-length queues, the selection of commands which result in errors seems a bit "random".

Your pde is first-order in time and third-order in space. To solve it you are (probably) going to need one initial condition, which you have (ie u(x,0) = 2 / cosh2(x) ) and three boundary conditions, which you don't have - something like u(0, t)=f(t), u(10,t)=g(t), D(u)(0,t)=h(t). I suggest you try to come up with the three relevant boundary conditions.

@Preben Alsholm 

I'm running Maple 2022 on Windows 7 - see attached

  restart:
  interface(version);

`Standard Worksheet Interface, Maple 2022.0, Windows 7, March 8 2022 Build ID 1599809`

(1)

#
# This works
#
  A:=DEQueue();
  empty(A);  
  push_back(A,2);
  

module DEQueue () local num, head, tail, storage, dsp; option object; end module

 

true

 

module DEQueue () local num, head, tail, storage, dsp; option object; end module

(2)

#
# Bit this results in an error!
#
  B:=DEQueue();
  empty(B);
  push_front(B,2);

module DEQueue () local num, head, tail, storage, dsp; option object; end module

 

true

 

Error, invalid return value from method moduledefinition: 'NULL'

 

 

Download DEQueue_Prob2.mw

@jjweimer 

You state

Maple -> PDF is a kludge whether on macOS or Windows, but the kludge on Windows works a bit better.

Actually, on Windows it works all the time, every time and has done for many years - not sure I'd call that a "kludge" The current (2022) Maple help for this exercise states

Windows and Linux
1. Open the worksheet that you want to export.
2. From the File menu, select Export As. The Export As dialog opens.
3. Select PDF as a file type.
4. Specify a path and folder for the file.
5. Enter a filename.
6. Click Save.

The same help for Maple 18 (ie released in ~2014) states

Windows, UNIX, and Linux
1. Open the worksheet that you want to export.
2. From the File menu, select Export As. The Export As dialog opens.
3. Select PDF as a file type.
4. Specify a path and folder for the file.
5. Enter a filename.
Click Save.

So what you refer to as a "kludge" has been working perfectly for at least the last seven years. I'd actually say a lot longer, because I have been using it for longer. Unfortunately I haven't got any Maple releases earlier than Maple 18 (circa 2014) loaded on my current machine, so I can't verify. As stated in my original response - I've never tried this on a MAC.

Can only suggest that you contact Maple Technical support, rather than the user community here

@michalkvasnicka 

  1. You can't install Updates version 1183 into Maple 2021.2, the most recent version which is supported in this Maple release is 1181
  2. Slightly puzzled as to why in Maple 2022, it appears to be looking in  the directory /home/kva/maple/toolbox/2021/Physics Updates/lib/ it should be looking in the directory /home/kva/maple/toolbox/2022/Physics Updates/lib/ it

What I'd do at this point

  1. Close all Maple sessions
  2. Go to the directory /home/kva/maple/toolbox/2021  and delete the Physics Updates folder
  3. Go to the directory /home/kva/maple/toolbox/2022  and delete the Physics Updates folder
  4. Start Maple 2022, and enter Physics:-Version(latest), close Maple
  5. Start Maple 2021 and enter Physics:-Version(1181), close Maple
  6. Start  both versions again and see what you have with Physics:-Version()

@ijuptilk 

  1. Firstly Acer's answer with Student[Calcculus1]:-Roots() was "better" than mine because it found more roots, although see item 2 below for some concerns that I have. My "installation" of the DirectSearch package, simply consists of having three files DirectSearch.hdb, DirectSearch.helpcand DirectSearch.mla in the directory (Windows machine) C:\Program Files\Maple 2022\lib. If you have these files in this directory, then close Maple and start again. I'm pretty sure that the 'lib' diirectory is only read on start-up. If you installed the DirectSearch() package somewhere else, then you are going to have to edit the 'libnam'e variable to include the location (see ?libname).
  2. Your specific roots problem is interesting, just because the expression is "evil". I was comparing the results I had obtained with Acer's results to see which roots DirectSeach had missed (Acer had 32 positive roots I had 21), just to see if I could figure out why. Moderately careful examination of the roots Acer provided, suggested that these were in fact all odd multiples of Pi/2. Looking at the second term (the fractional one) in your expression, both numerator and denominator contain tan(q) which is infinite for all odd multiples of Pi/2 - so this ought to lead to infinity divided by infinity, which I would consider "undefined". If one actually evaluates your expression at an odd multiple of Pi/2, then Maple generates a numeric exception error.; So how could these values possibly be "roots". In fact my original assumption was incorrect, the roots are not quite at odd multiples of Pi/2, but are very close to these values, and get closer as the value of the root increases (see the attached). I can only imagine that the proximity of the actual roots to values where the function is undefined is causing significant problems for any rootfinding algorithm - maybe this is why the different approaches have such different levesl of success?? Using RootFinding:-NextZero() with a bit of guidance, I could produce 6 or 7 roots, fsolve() with the 'avoid' option only ever produces 2, DirectSearch() produces 21 positive roots and

** Rant mode on**
I still think that Maple should do something more intelligent than throw a numeric exception when given tan(Pi/2). I vaguely remember submitting an SCR about this many years ago. This behaviour has occured in the last 10 major Maple releases - all the way back to Maple 18 (and maybe before)
**Rant Mode off**

restart:

with(Student:-Calculus1):

a[1] := .1093: k[3] := 7.5*10^(-12): k[2] := 3.8*10^(-12):
d := 0.2e-3: eta[1] := 0.240e-1: alpha[2] := -.1104:
alpha[3] := -0.1104e-2: eta[2] := .1361: xi := 1.219*10^(-6):
alpha := 1-alpha[3]^2/(a[1]*eta[1]):theta[0] := 0.5e-1:
Hc := (Pi/d)*sqrt(k[2]/xi): H := 5.5*Hc:
lambda := a[1]/(xi*H^2):

               
sols := Roots((H/(Hc))^2 -(4*q^2)/Pi^2*((tan(q)- q/(1-alpha))/(tan(q)-q)),
              q=0..100, numeric, maxsols=200);

[1.579423766, 4.712622100, 7.853994064, 10.99555811, 14.13714635, 17.27873942, 20.42033357, 23.56192783, 26.70352198, 29.84511595, 32.98670974, 36.12830339, 39.26989691, 42.41149032, 45.55308364, 48.69467689, 51.83627007, 54.97786320, 58.11945627, 61.26104931, 64.40264231, 67.54423529, 70.68582823, 73.82742115, 76.96901405, 80.11060693, 83.25219980, 86.39379265, 89.53538549, 92.67697831, 95.81857112, 98.96016393]

(1)

#
# Check distance between roots and odd multiples of Pi/2
#
  seq( sols[j]-(2*j-1)*Pi/2, j=1..numelems(sols));

0.8627439e-2, 0.233119e-3, 0.12429e-4, -0.1618e-4, -0.2059e-4, -0.2018e-4, -0.1868e-4, -0.1707e-4, -0.1558e-4, -0.1426e-4, -0.1313e-4, -0.1213e-4, -0.1127e-4, -0.1051e-4, -0.984e-5, -0.925e-5, -0.872e-5, -0.824e-5, -0.783e-5, -0.744e-5, -0.710e-5, -0.677e-5, -0.649e-5, -0.622e-5, -0.597e-5, -0.575e-5, -0.553e-5, -0.533e-5, -0.515e-5, -0.498e-5, -0.483e-5, -0.467e-5

(2)

 

 

Download roots3.mw

@ijuptilk 

it depends on how you want the answers to be stored for subsequent access - eg Matrix, Array, table, list, Vector whatever, which in turn, probably deends on what you are going to do wqith these valkue. The atached compues the function you specify at each of the root values and stores the result in an nx1 matrix. You want it store this in any other kind of "data container", thsat would be pretty trivial to do

  restart:
  with(RootFinding):
  a[1]:= .1093: k[3]:= 7.5*10^(-12): k[2]:= 3.8*10^(-12):
  d:= 0.2e-3: eta[1]:= 0.240e-1: alpha[2]:= -.1104:
  alpha[3]:= -0.1104e-2: eta[2]:= .1361: xi:= 1.219*10^(-6):
  alpha:= 1-alpha[3]^2/(a[1]*eta[1]): theta[0]:= 0.5e-1:
  Hc:= (Pi/d)*sqrt(k[2]/xi):

  H:= 5.5*Hc:
  lambda:= a[1]/(xi*H^2):
  f:= unapply((H/Hc)^2-4*q^2*(tan(q)-q/(1-alpha))/(Pi^2*(tan(q)-q)),q);

proc (q) options operator, arrow; 30.24999998-4*q^2*(tan(q)-2152.252494*q)/(Pi^2*(tan(q)-q)) end proc

(1)

  interface(rtablesize=50):
  ans:=DirectSearch:-SolveEquations(f,AllSolutions=true);
#
# Extract the actual roots and sort them
#
  ansSort:=sort(entries~(ans[..,3], 'nolist'));
#
# Compute some function with these roots
#
  p:=Matrix( 1..numelems(ansSort),1, i->alpha*lambda/(1-(4*ansSort[i])^2*(Hc/H)^2/Pi^2) );
  interface(rtablesize=10):

Matrix(47, 4, {(1, 1) = 0.1609764636e-15, (1, 2) = Vector(1, {(1) = 0.1268765004e-7}), (1, 3) = Vector(1, {(1) = -1.5794237664}), (1, 4) = 21, (2, 1) = 0.1590985127e-14, (2, 2) = Vector(1, {(1) = -0.3988715491e-7}), (2, 3) = Vector(1, {(1) = 1.5794237664}), (2, 4) = 42, (3, 1) = 0.1360469941e-11, (3, 2) = Vector(1, {(1) = -0.11664e-5}), (3, 3) = Vector(1, {(1) = 42.4114903226}), (3, 4) = 59, (4, 1) = 0.1757286539e-11, (4, 2) = Vector(1, {(1) = 0.13256e-5}), (4, 3) = Vector(1, {(1) = -89.5353854863}), (4, 4) = 67, (5, 1) = 0.8651284813e-11, (5, 2) = Vector(1, {(1) = -0.29413e-5}), (5, 3) = Vector(1, {(1) = 17.2787394247}), (5, 4) = 47, (6, 1) = 0.1688036992e-10, (6, 2) = Vector(1, {(1) = 0.41086e-5}), (6, 3) = Vector(1, {(1) = -86.3937926494}), (6, 4) = 55, (7, 1) = 0.5109618415e-10, (7, 2) = Vector(1, {(1) = 0.71482e-5}), (7, 3) = Vector(1, {(1) = 80.1106069341}), (7, 4) = 70, (8, 1) = 0.5109618415e-10, (8, 2) = Vector(1, {(1) = 0.71482e-5}), (8, 3) = Vector(1, {(1) = -80.1106069341}), (8, 4) = 69, (9, 1) = 0.6767702417e-10, (9, 2) = Vector(1, {(1) = 0.82266e-5}), (9, 3) = Vector(1, {(1) = -73.8274211520}), (9, 4) = 69, (10, 1) = 0.1035917565e-9, (10, 2) = Vector(1, {(1) = 0.101780e-4}), (10, 3) = Vector(1, {(1) = 29.8451159450}), (10, 4) = 47, (11, 1) = 0.1038611085e-9, (11, 2) = Vector(1, {(1) = 0.101912e-4}), (11, 3) = Vector(1, {(1) = -4.7126220995}), (11, 4) = 27, (12, 1) = 0.1072475664e-9, (12, 2) = Vector(1, {(1) = 0.103560e-4}), (12, 3) = Vector(1, {(1) = 26.7035219764}), (12, 4) = 51, (13, 1) = 0.3458136494e-9, (13, 2) = Vector(1, {(1) = 0.185961e-4}), (13, 3) = Vector(1, {(1) = -36.1283033908}), (13, 4) = 48, (14, 1) = 0.3584630775e-9, (14, 2) = Vector(1, {(1) = 0.189331e-4}), (14, 3) = Vector(1, {(1) = -17.2787394247}), (14, 4) = 40, (15, 1) = 0.3658928321e-9, (15, 2) = Vector(1, {(1) = 0.191283e-4}), (15, 3) = Vector(1, {(1) = 89.5353854863}), (15, 4) = 56, (16, 1) = 0.6622298869e-9, (16, 2) = Vector(1, {(1) = -0.257338e-4}), (16, 3) = Vector(1, {(1) = -64.4026423139}), (16, 4) = 66, (17, 1) = 0.6703036941e-9, (17, 2) = Vector(1, {(1) = 0.258902e-4}), (17, 3) = Vector(1, {(1) = -23.5619278325}), (17, 4) = 46, (18, 1) = 0.6900373055e-9, (18, 2) = Vector(1, {(1) = -0.262686e-4}), (18, 3) = Vector(1, {(1) = -32.9867097431}), (18, 4) = 52, (19, 1) = 0.7288453593e-9, (19, 2) = Vector(1, {(1) = 0.269971e-4}), (19, 3) = Vector(1, {(1) = 4.7126220993}), (19, 4) = 27, (20, 1) = 0.7667552900e-9, (20, 2) = Vector(1, {(1) = 0.276903e-4}), (20, 3) = Vector(1, {(1) = -7.8539940639}), (20, 4) = 32, (21, 1) = 0.7701963329e-9, (21, 2) = Vector(1, {(1) = -0.277524e-4}), (21, 3) = Vector(1, {(1) = 76.9690140524}), (21, 4) = 71, (22, 1) = 0.1018828140e-8, (22, 2) = Vector(1, {(1) = -0.319191e-4}), (22, 3) = Vector(1, {(1) = 7.8539940641}), (22, 4) = 43, (23, 1) = 0.1054823577e-8, (23, 2) = Vector(1, {(1) = -0.324780e-4}), (23, 3) = Vector(1, {(1) = -45.5530836440}), (23, 4) = 66, (24, 1) = 0.1086951117e-8, (24, 2) = Vector(1, {(1) = -0.329689e-4}), (24, 3) = Vector(1, {(1) = 54.9778631952}), (24, 4) = 60, (25, 1) = 0.1789054721e-8, (25, 2) = Vector(1, {(1) = -0.422972e-4}), (25, 3) = Vector(1, {(1) = 64.4026423139}), (25, 4) = 56, (26, 1) = 0.2414497685e-8, (26, 2) = Vector(1, {(1) = 0.491375e-4}), (26, 3) = Vector(1, {(1) = 51.8362700697}), (26, 4) = 64, (27, 1) = 0.2793647501e-8, (27, 2) = Vector(1, {(1) = 0.528550e-4}), (27, 3) = Vector(1, {(1) = 58.1194562736}), (27, 4) = 52, (28, 1) = 0.3021817696e-8, (28, 2) = Vector(1, {(1) = -0.549711e-4}), (28, 3) = Vector(1, {(1) = 61.2610493114}), (28, 4) = 46, (29, 1) = 0.3039225144e-8, (29, 2) = Vector(1, {(1) = -0.551292e-4}), (29, 3) = Vector(1, {(1) = 23.5619278325}), (29, 4) = 52, (30, 1) = 0.3453070599e-8, (30, 2) = Vector(1, {(1) = -0.587628e-4}), (30, 3) = Vector(1, {(1) = 48.6946768892}), (30, 4) = 62, (31, 1) = 0.3582400012e-8, (31, 2) = Vector(1, {(1) = 0.598532e-4}), (31, 3) = Vector(1, {(1) = -42.4114903226}), (31, 4) = 67, (32, 1) = 0.3734419284e-8, (32, 2) = Vector(1, {(1) = 0.611099e-4}), (32, 3) = Vector(1, {(1) = -29.8451159450}), (32, 4) = 49, (33, 1) = 0.5343604629e-8, (33, 2) = Vector(1, {(1) = 0.731000e-4}), (33, 3) = Vector(1, {(1) = 73.8274211520}), (33, 4) = 62, (34, 1) = 0.5902038576e-8, (34, 2) = Vector(1, {(1) = 0.768247e-4}), (34, 3) = Vector(1, {(1) = -26.7035219764}), (34, 4) = 54, (35, 1) = 0.6032306911e-8, (35, 2) = Vector(1, {(1) = 0.776679e-4}), (35, 3) = Vector(1, {(1) = -67.5442352858}), (35, 4) = 61, (36, 1) = 0.7371464447e-8, (36, 2) = Vector(1, {(1) = 0.858572e-4}), (36, 3) = Vector(1, {(1) = -48.6946768892}), (36, 4) = 63, (37, 1) = 0.8290250265e-8, (37, 2) = Vector(1, {(1) = -0.910508e-4}), (37, 3) = Vector(1, {(1) = 36.1283033908}), (37, 4) = 58, (38, 1) = 0.8477611289e-8, (38, 2) = Vector(1, {(1) = 0.920739e-4}), (38, 3) = Vector(1, {(1) = 86.3937926494}), (38, 4) = 64, (39, 1) = 0.1071247208e-7, (39, 2) = Vector(1, {(1) = 0.1035011e-3}), (39, 3) = Vector(1, {(1) = 45.5530836440}), (39, 4) = 61, (40, 1) = 0.1095995057e-7, (40, 2) = Vector(1, {(1) = 0.1046898e-3}), (40, 3) = Vector(1, {(1) = 39.2698969106}), (40, 4) = 47, (41, 1) = 0.1283914116e-7, (41, 2) = Vector(1, {(1) = 0.1133099e-3}), (41, 3) = Vector(1, {(1) = -70.6858282308}), (41, 4) = 57, (42, 1) = 0.1286114994e-7, (42, 2) = Vector(1, {(1) = -0.1134070e-3}), (42, 3) = Vector(1, {(1) = 67.5442352858}), (42, 4) = 67, (43, 1) = 0.1480693513e-7, (43, 2) = Vector(1, {(1) = 0.1216838e-3}), (43, 3) = Vector(1, {(1) = -51.8362700697}), (43, 4) = 61, (44, 1) = 0.1863711818e-7, (44, 2) = Vector(1, {(1) = -0.1365178e-3}), (44, 3) = Vector(1, {(1) = -20.4203335662}), (44, 4) = 52, (45, 1) = 0.1939169567e-7, (45, 2) = Vector(1, {(1) = 0.1392541e-3}), (45, 3) = Vector(1, {(1) = -98.9601639287}), (45, 4) = 55, (46, 1) = 0.2490496854e-7, (46, 2) = Vector(1, {(1) = -0.1578131e-3}), (46, 3) = Vector(1, {(1) = -76.9690140524}), (46, 4) = 62, (47, 1) = 0.3092106349e-7, (47, 2) = Vector(1, {(1) = 0.1758439e-3}), (47, 3) = Vector(1, {(1) = -92.6769783110}), (47, 4) = 54})

 

Vector(47, {(1) = -98.9601639287, (2) = -92.6769783110, (3) = -89.5353854863, (4) = -86.3937926494, (5) = -80.1106069341, (6) = -76.9690140524, (7) = -73.8274211520, (8) = -70.6858282308, (9) = -67.5442352858, (10) = -64.4026423139, (11) = -51.8362700697, (12) = -48.6946768892, (13) = -45.5530836440, (14) = -42.4114903226, (15) = -36.1283033908, (16) = -32.9867097431, (17) = -29.8451159450, (18) = -26.7035219764, (19) = -23.5619278325, (20) = -20.4203335662, (21) = -17.2787394247, (22) = -7.8539940639, (23) = -4.7126220995, (24) = -1.5794237664, (25) = 1.5794237664, (26) = 4.7126220993, (27) = 7.8539940641, (28) = 17.2787394247, (29) = 23.5619278325, (30) = 26.7035219764, (31) = 29.8451159450, (32) = 36.1283033908, (33) = 39.2698969106, (34) = 42.4114903226, (35) = 45.5530836440, (36) = 48.6946768892, (37) = 51.8362700697, (38) = 54.9778631952, (39) = 58.1194562736, (40) = 61.2610493114, (41) = 64.4026423139, (42) = 67.5442352858, (43) = 73.8274211520, (44) = 76.9690140524, (45) = 80.1106069341, (46) = 86.3937926494, (47) = 89.5353854863})

 

Matrix(%id = 36893488148166952948)

(2)

 

 

Download roots2.mw

@nahid200

Why am I fixing this crap code?

  restart:  
  Wpath:= proc(steps,t)
               local walk, i, N, ww, j:
               N := nops(steps):
               walk [0]:= 0:  
               for i from 0 to N-1 do  
                   walk[i+1]:= walk[i]+ steps[i+1]*sqrt(t/N)*stats[random,normald[0,1]](1):
               od:
               plot([seq( [j, walk[j]], j=0..N)]);
          end proc:

  Wpath([seq(j, j=1..10^3)], 1);

 

 

 

Download RW.mw

@AHSAN 

for wi=0.1, k=0.1, x=0.0, Q=0.5516, lambda=-0.0041. I get the attached. Fit isn't quite as good but not bad.

odeplots7.mw

First 20 21 22 23 24 25 26 Last Page 22 of 207