Advanced Search

MaplePrimes Questions

fix the collect ...

Asked by: bashar27 55
user-error recursive-assignment duplicate-question
September 11 2019
0  2

restart;
T := mu+lambda*H(xi)+(v-1)*H(xi)^2;
                                                2
               mu + lambda H(xi) + (v - 1) H(xi) 
u[0] := a[0]+a[1]*(d+H(xi))+a[2]/(d+H(xi))+a[3]*(d+H(xi))^2+a[4]/(d+H(xi))^2;
                                a[2]                      2
    a[0] + a[1] (d + H(xi)) + --------- + a[3] (d + H(xi)) 
                              d + H(xi)                    

             a[4]    
       + ------------
                    2
         (d + H(xi)) 
diff(u[0], xi);
                            / d        \
                       a[2] |---- H(xi)|
        / d        \        \ dxi      /
   a[1] |---- H(xi)| - -----------------
        \ dxi      /                2   
                         (d + H(xi))    

                                                 / d        \
                                          2 a[4] |---- H(xi)|
                           / d        \          \ dxi      /
      + 2 a[3] (d + H(xi)) |---- H(xi)| - -------------------
                           \ dxi      /                 3    
                                             (d + H(xi))     
collect(%, diff(H(xi), xi));
/           a[2]                               2 a[4]   \ / d       
|a[1] - ------------ + 2 a[3] (d + H(xi)) - ------------| |---- H(xi
|                  2                                   3| \ dxi     
\       (d + H(xi))                         (d + H(xi)) /           

   \
  )|
   /
d[1] := (a[1]-a[2]/(d+H(xi))^2+2*a[3]*(d+H(xi))-2*a[4]/(d+H(xi))^3)*T;
 /           a[2]                               2 a[4]   \ /  
 |a[1] - ------------ + 2 a[3] (d + H(xi)) - ------------| \mu
 |                  2                                   3|    
 \       (d + H(xi))                         (d + H(xi)) /    

                                  2\
    + lambda H(xi) + (v - 1) H(xi) /
diff(d[1], xi);
/       / d        \                                / d        \\ 
|2 a[2] |---- H(xi)|                         6 a[4] |---- H(xi)|| 
|       \ dxi      /          / d        \          \ dxi      /| 
|------------------- + 2 a[3] |---- H(xi)| + -------------------| 
|              3              \ dxi      /                 4    | 
\   (d + H(xi))                                 (d + H(xi))     / 

  /                                 2\   /           a[2]    
  \mu + lambda H(xi) + (v - 1) H(xi) / + |a[1] - ------------
                                         |                  2
                                         \       (d + H(xi)) 

                             2 a[4]   \ /       / d        \
   + 2 a[3] (d + H(xi)) - ------------| |lambda |---- H(xi)|
                                     3| \       \ dxi      /
                          (d + H(xi)) /                     

                     / d        \\
   + 2 (v - 1) H(xi) |---- H(xi)||
                     \ dxi      //
collect(%, diff(H(xi), xi));
//   2 a[2]                  6 a[4]   \ /                 
||------------ + 2 a[3] + ------------| \mu + lambda H(xi)
||           3                       4|                   
\\(d + H(xi))             (d + H(xi)) /                   

                  2\   /           a[2]                         
   + (v - 1) H(xi) / + |a[1] - ------------ + 2 a[3] (d + H(xi))
                       |                  2                     
                       \       (d + H(xi))                      

        2 a[4]   \                           \ / d        \
   - ------------| (lambda + 2 (v - 1) H(xi))| |---- H(xi)|
                3|                           | \ dxi      /
     (d + H(xi)) /                           /             
d[2] := ((2*a[2]/(d+H(xi))^3+2*a[3]+6*a[4]/(d+H(xi))^4)*(mu+lambda*H(xi)+(v-1)*H(xi)^2)+(a[1]-a[2]/(d+H(xi))^2+2*a[3]*(d+H(xi))-2*a[4]/(d+H(xi))^3)*(lambda+(2*(v-1))*H(xi)))*T;
//   2 a[2]                  6 a[4]   \ /                 
||------------ + 2 a[3] + ------------| \mu + lambda H(xi)
||           3                       4|                   
\\(d + H(xi))             (d + H(xi)) /                   

                  2\   /           a[2]                         
   + (v - 1) H(xi) / + |a[1] - ------------ + 2 a[3] (d + H(xi))
                       |                  2                     
                       \       (d + H(xi))                      

        2 a[4]   \                           \ /                 
   - ------------| (lambda + 2 (v - 1) H(xi))| \mu + lambda H(xi)
                3|                           |                   
     (d + H(xi)) /                           /                   

                  2\
   + (v - 1) H(xi) /

eq := (2*k*k)*w*beta*d[2]-(2*alpha*k*k)*d[1]-2*w*u[0]+k*u[0]*u[0];
   2        //   2 a[2]                  6 a[4]   \ /  
2 k  w beta ||------------ + 2 a[3] + ------------| \mu
            ||           3                       4|    
            \\(d + H(xi))             (d + H(xi)) /    

                                 2\   /           a[2]    
   + lambda H(xi) + (v - 1) H(xi) / + |a[1] - ------------
                                      |                  2
                                      \       (d + H(xi)) 

                             2 a[4]   \                          
   + 2 a[3] (d + H(xi)) - ------------| (lambda + 2 (v - 1) H(xi)
                                     3|                          
                          (d + H(xi)) /                          

   \ /                                 2\            2 /    
  )| \mu + lambda H(xi) + (v - 1) H(xi) / - 2 alpha k  |a[1]
   |                                                   |    
   /                                                   \    

         a[2]                               2 a[4]   \ /  
   - ------------ + 2 a[3] (d + H(xi)) - ------------| \mu
                2                                   3|    
     (d + H(xi))                         (d + H(xi)) /    

                                 2\       /    
   + lambda H(xi) + (v - 1) H(xi) / - 2 w |a[0]
                                          |    
                                          \    

                          a[2]                      2
   + a[1] (d + H(xi)) + --------- + a[3] (d + H(xi)) 
                        d + H(xi)                    

         a[4]    \     /                            a[2]   
   + ------------| + k |a[0] + a[1] (d + H(xi)) + ---------
                2|     |                          d + H(xi)
     (d + H(xi)) /     \                                   

                     2       a[4]    \  
   + a[3] (d + H(xi))  + ------------|^2
                                    2|  
                         (d + H(xi)) /  
value(%);
   2        //   2 a[2]                  6 a[4]   \ /  
2 k  w beta ||------------ + 2 a[3] + ------------| \mu
            ||           3                       4|    
            \\(d + H(xi))             (d + H(xi)) /    

                                 2\   /           a[2]    
   + lambda H(xi) + (v - 1) H(xi) / + |a[1] - ------------
                                      |                  2
                                      \       (d + H(xi)) 

                             2 a[4]   \                          
   + 2 a[3] (d + H(xi)) - ------------| (lambda + 2 (v - 1) H(xi)
                                     3|                          
                          (d + H(xi)) /                          

   \ /                                 2\            2 /    
  )| \mu + lambda H(xi) + (v - 1) H(xi) / - 2 alpha k  |a[1]
   |                                                   |    
   /                                                   \    

         a[2]                               2 a[4]   \ /  
   - ------------ + 2 a[3] (d + H(xi)) - ------------| \mu
                2                                   3|    
     (d + H(xi))                         (d + H(xi)) /    

                                 2\       /    
   + lambda H(xi) + (v - 1) H(xi) / - 2 w |a[0]
                                          |    
                                          \    

                          a[2]                      2
   + a[1] (d + H(xi)) + --------- + a[3] (d + H(xi)) 
                        d + H(xi)                    

         a[4]    \     /                            a[2]   
   + ------------| + k |a[0] + a[1] (d + H(xi)) + ---------
                2|     |                          d + H(xi)
     (d + H(xi)) /     \                                   

                     2       a[4]    \  
   + a[3] (d + H(xi))  + ------------|^2
                                    2|  
                         (d + H(xi)) /  
expr := simplify(%);
Error, (in simplify) too many levels of recursion
temp := algsubs(d+H(xi) = freeze(d+H(xi)), numer(expr));
                              expr
thaw(collect(temp, freeze(d+H(xi)))/denom(expr));
                              expr
collect(%, H(xi));
 

Why certain delimeters act in a strange way...

Asked by: Adam Ledger 360
string sprintf escape
September 11 2019
0  2

Can someone please explain to me why this occurs:

 


 

> 

with(StringTools):
> 

Join(["H:\\USB 1 BACKUP\\ESD-USB\\", "Chemical Engineering"])

"H:\USB 1 BACKUP\ESD-USB\ Chemical Engineering"

 (1)
> 

convert("H:\\USB 1 BACKUP\\ESD-USB\\ Chemical Engineering", 'symbol')

`H:\USB 1 BACKUP\ESD-USB\ Chemical Engineering`

 (2)
> 

convert('`H:\USB 1 BACKUP\ESD-USB\ Chemical Engineering`', 'string')

"H:USB 1 BACKUPESD-USBChemical Engineering"

 (3)
> 

``


 

Download this_makes_me_grumpy.mw

Solving inequalities with conditions...

Asked by: vanzzy 35
solve inequality
September 11 2019
0  2

Hello everyone,

I'm struggling to solve inequalities with conditions.

I have this inequality with 4 variables, which I have some conditions. However, I can't implement this conditions to the inequality and solve using the 'solve' command.

Can anybody help me?

inequality.mw

What does VectorCalculus:-Curvature represent?...

Asked by: sand15 812
curvature vectorcalculus
September 11 2019
0  4

Hi,

I'm surprised by the result of the procedure VectorCalculus:-Curvature which is always a positive scalar quantity:
For instance
c := VectorCalculus:-Curvature(<x, sin(x)>, x):
plot(c, x=0..2*Pi)  # c >=0 for all x in [0, 2*Pi]

In the help pages it's written that the (signed) curvature for a function y(x) is y''/(1+y' 2)(3/2).

y := sin(x):
c := diff(y, x$2) / (1+diff(y,x)^2)^(3/2):
plot(c,  x=0..2*Pi)  # c < 0 if x in (0, Pi)  and  c > 0 if x in (Pi, 2*Pi)

Could you please help me to understand this?

Thank in advance

When can we expect enhancement of Heun functions?...

Asked by: J F Ogilvie 258
special-function differential-equations heun
September 10 2019
1  6

Some years ago it was promised that expansion of capabilities of Heun functions was imminent, but nothing has appeared.  Other functions long overdue for inclusion as special functions in Maple are the Lame functions, which arise as special cases of Heun's differential equation and therefore of Heun functions.  Lame's differential equation appears in Abramowitz and Stegun, but has long been neglected in Maple.  These spectial functions are much more generally useful to users of Maple than, for instance, esoteric parts of the physics package. 

Starting a new java session...

Asked by: denitsastaicova 80
gui java
September 10 2019
3  2

I remember there was a command that used to start an independent java session each time. I thought it's "xmaple -singlemode" but it's not. So could anybody please remind me what was the command. I just killed a second java session after losing kernel connection and java running wild. I'd be very happy to do my calculation in another java session so that if anything goes bad II won't have to open all my 5 tabs again.

How to convert a file from mathematica notebook t...

Asked by: Jayaprakash J 25
mathematica import duplicate-question
September 10 2019
1  7

tittle.6.mw 
 

> 

"Manipulate[Module[{sol,p1,p2,sys,tbleq,eq,init,sol1,solDTM,xDTM}, Np=10;H=1; sol=NDSolve[{alpha D[u[x,t],{x,2}]== D[u[x,t],{t,1}],u[-1,t]== 1,u[1,t]== 1,u[x,0]== Exp[-1000 (x+1)^2 (1-x)^2]},u[x,t],{x,-1,1},{t,0,1}]; p1=Plot[u[x,t]/.sol/.t-> ta,{x,-1,1},Frame-> True,PlotStyle->{Red,Thick},GridLines->Automatic,PlotRange->{0,1}]; eq[k_]:=alpha (k+2) (k+1) F[k+2][t]- F[k]'[t]==0; tbleq=Table[eq[i],{i,0,Np}]; init=Table[F[k][0]== If[k==0,0,0],{k,0,Np}]; sol1=Solve[{Sum[(1/H)^k F[k][t],{k,0,Np+2}]==1,Sum[(-1/H)^k F[k][t],{k,0,Np+2}]==1},{F[Np+1][t],F[Np+2][t]}]; sys=Join[tbleq/.sol1,init]//Flatten; solDTM=NDSolve[sys//Flatten,Table[F[k],{k,0,Np}],{t,0,0.5}]; xDTM[x_]:= Sum[((x/H)^k) F[k][t]/.sol1/.solDTM,{k,0,Np+2}]; p2=Plot[xDTM[x]/.t-> ta,{x,-1,1},Frame->True,GridLines->Automatic,PlotStyle->{Dashed,Thick,Blue},PlotRange->{0,1}]; Show[p1,p2,ImageSize->{580,400}]],{{alpha,1,"thermal diffusivity"},0.1,2,0.1,Appearance-> "Labeled"},{{ta,0.075,"time"},0.01,0.1,0.005,Appearance-> "Labeled"}] "


 

Download tittle.6.mw

sum, diff and subs...

Asked by: zhuxian 60
sum diff subs
September 09 2019
0  3

I am sorry. There is no 'tag' button on my computer, so I can't ask question here. I don't know the reason.

So I ask here. Now, my real question is:

sum(subs(x = 0, diff(x^2, x $ k)), k = 1 .. 2);

                               0

subs(x = 0, diff(x^2, x $ 2));
                               2


I think the first result should be 2, but it gives me 0.

Thanks in advance.

 

Problem in finding the solution of linear algebrai...

Asked by: Muhammad Usman 240
fsolve
September 09 2019
0  2

Dear Users!

Hope you would be fine with everything. I want to find the solution of linear algebric equations but fsolve command not working please see and fix this problem. I shall be very thankful.

C[0] := 3.19153824321146142351956847947*tau[1]-19.1492294592687685411174108768*tau[2]+111.703838512401149823184896781*tau[3]+3.19153824321146142351956847947*tau[4]-44.6815354049604599292739587124*tau[5]+622.349957426234977586315853494*tau[6];
C[1] := 51.0646118913833827763130956714*tau[2]-612.775342696600593315757148056*tau[3]+51.0646118913833827763130956714*tau[5]-1429.80913295873471773676667880*tau[6];
C[2] := -1.06073680388443795908856507616+3.19153824321146142351956847947*tau[1]+53.1609155734306093706448370717*tau[2]+1672.89412862088744108725223170*tau[3]+3.19153824321146142351956847947*tau[4]+27.6286096277389179824882892361*tau[5]+1026.57792701153122226218722129*tau[6];
C[3] := -1.08847004231036963538035920033+3.19153824321146142351956847947*tau[1]+62.6399144226357196540662623767*tau[2]+2040.52109049201342887896297462*tau[3]+3.19153824321146142351956847947*tau[4]+37.1076084769440282659097145411*tau[5]+1242.54090729537544551915515930*tau[6];
C[4] := -1.05523181556926815105314303389+3.19153824321146142351956847947*tau[1]+72.7671212023804312453829273862*tau[2]+2472.93216226733267613216245895*tau[3]+3.19153824321146142351956847947*tau[4]+47.2348152566887398572263795506*tau[5]+1512.91667059477930731128800348*tau[6];
C[5] := -.922876006485286011069063957991+3.19153824321146142351956847947*tau[1]+82.9822841707707093164204255644*tau[2]+2971.36790137532483139495115633*tau[3]+3.19153824321146142351956847947*tau[4]+57.4499782250790179282638777288*tau[5]+1847.90980220852701343747673000*tau[6];

fsolve({seq(`$`(C[l1], l1 = 0 .. 5))});

Special request to:
@acer @Carl Love @Kitonum @Preben Alsholm

Terrible computations involving eulers number...

Asked by: mschneider 60
numeric evalf user-error
September 09 2019
0  9

Please help...

answer := 5*exp(7.5*t);

response := subs(a = exp(1), 5*a^(7.5*t));

a := evalf(subs(t = Pi, answer));

b := evalf(subs(t = Pi, response));

evalf(a-b);

 

Any explanations as to why the last line is not zero?  any workarounds?

fix simplify problem...

Asked by: bashar27 55
September 09 2019
0  2

restart;
lambda := 1;
                               1
mu := 1;
                               1
v := 2;
                               2
r := lambda*(v-1);
                               1
g := mu*(v-1);
                               1
a[2] := 0;
                               0

omega := -(1/2)*alpha*l[1]*l[2]*lambda^2+2*alpha*l[1]*l[2]*mu*v-2*alpha*l[1]*l[2]*mu-alpha*h[1]*h[2];
              3                                  
              - alpha l[1] l[2] - alpha h[1] h[2]
              2                                  
a[0] := -(1/2)*(2*d*v-2*d-lambda)*alpha*l[1]*l[2]/(h[1]*beta*(sqrt(h[1]*beta*alpha*l[1]*l[2])/(h[1]*beta)));
                    (2 d - 1) alpha l[1] l[2]     
              - ----------------------------------
                                             (1/2)
                2 (h[1] beta alpha l[1] l[2])     
a[1] := sqrt(h[1]*beta*alpha*l[1]*l[2])*(v-1)/(beta*h[1]);
                                           (1/2)
                (h[1] beta alpha l[1] l[2])     
                --------------------------------
                           beta h[1]            

Omega := lambda^2-4*mu*v+4*mu;
                               -3
H := (-lambda+sqrt(-Omega)*{tan(sqrt(-Omega)*xi)+sec(sqrt(-Omega)*xi)})/(2*(v-1));
         1   1  (1/2)  /   / (1/2)   \      / (1/2)   \\ 
       - - + - 3      { tan\3      xi/ + sec\3      xi/ }
         2   2         \                               / 

u := a[0]+a[1]*(d+H)+a[2]/(d+H);
            (2 d - 1) alpha l[1] l[2]            1     /
      - ---------------------------------- + --------- |
                                     (1/2)   beta h[1] \
        2 (h[1] beta alpha l[1] l[2])                   

                                   (1/2) /    1
        (h[1] beta alpha l[1] l[2])      |d - -
                                         \    2

           1  (1/2)  /   / (1/2)   \      / (1/2)   \\ \\
         + - 3      { tan\3      xi/ + sec\3      xi/ }||
           2         \                               / //
f := diff(u, xi);
               /                                         // 
        1      |                           (1/2)  (1/2) { | 
   ----------- \(h[1] beta alpha l[1] l[2])      3       \\1
   2 beta h[1]                                              

                      2\       
           / (1/2)   \ |  (1/2)
      + tan\3      xi/ / 3     

                                            \ \
           / (1/2)   \    / (1/2)   \  (1/2) }|
      + sec\3      xi/ tan\3      xi/ 3     / /
S := diff(f, xi);
             /                                         /      
      1      |                           (1/2)  (1/2) {      /
 ----------- \(h[1] beta alpha l[1] l[2])      3       \6 tan\
 2 beta h[1]                                                  

              /                  2\
    (1/2)   \ |       / (1/2)   \ |
   3      xi/ \1 + tan\3      xi/ /

                                     2
           / (1/2)   \    / (1/2)   \ 
    + 3 sec\3      xi/ tan\3      xi/ 

                       /                  2\\ \
           / (1/2)   \ |       / (1/2)   \ | }|
    + 3 sec\3      xi/ \1 + tan\3      xi/ // /

eq := -(alpha*h[1]*h[2]+omega)*u-(2*beta*h[1]*u*u)*u+alpha*l[1]*l[2]*S;
  3                 /      (2 d - 1) alpha l[1] l[2]        
- - alpha l[1] l[2] |- ---------------------------------- + 
  2                 |                               (1/2)   
                    \  2 (h[1] beta alpha l[1] l[2])        

      1     /                           (1/2) /    1
  --------- |(h[1] beta alpha l[1] l[2])      |d - -
  beta h[1] \                                 \    2

     1  (1/2)  /   / (1/2)   \      / (1/2)   \\ \\\             
   + - 3      { tan\3      xi/ + sec\3      xi/ }||| - 2 beta h[1
     2         \                               / //|             
                                                   /             

    /      (2 d - 1) alpha l[1] l[2]            1     /
  ] |- ---------------------------------- + --------- |
    |                               (1/2)   beta h[1] \
    \  2 (h[1] beta alpha l[1] l[2])                   

                             (1/2) /    1
  (h[1] beta alpha l[1] l[2])      |d - -
                                   \    2

     1  (1/2)  /   / (1/2)   \      / (1/2)   \\ \\\     
   + - 3      { tan\3      xi/ + sec\3      xi/ }|||^3 + 
     2         \                               / //|     
                                                   /     

              /                                                 
       1      |                                           (1/2) 
  ----------- \alpha l[1] l[2] (h[1] beta alpha l[1] l[2])      
  2 beta h[1]                                                   

          /                 /                  2\
   (1/2) {      / (1/2)   \ |       / (1/2)   \ |
  3       \6 tan\3      xi/ \1 + tan\3      xi/ /

                                    2
          / (1/2)   \    / (1/2)   \ 
   + 3 sec\3      xi/ tan\3      xi/ 

                      /                  2\\ \
          / (1/2)   \ |       / (1/2)   \ | }|
   + 3 sec\3      xi/ \1 + tan\3      xi/ // /
value(%);
  3                 /      (2 d - 1) alpha l[1] l[2]        
- - alpha l[1] l[2] |- ---------------------------------- + 
  2                 |                               (1/2)   
                    \  2 (h[1] beta alpha l[1] l[2])        

      1     /                           (1/2) /    1
  --------- |(h[1] beta alpha l[1] l[2])      |d - -
  beta h[1] \                                 \    2

     1  (1/2)  /   / (1/2)   \      / (1/2)   \\ \\\             
   + - 3      { tan\3      xi/ + sec\3      xi/ }||| - 2 beta h[1
     2         \                               / //|             
                                                   /             

    /      (2 d - 1) alpha l[1] l[2]            1     /
  ] |- ---------------------------------- + --------- |
    |                               (1/2)   beta h[1] \
    \  2 (h[1] beta alpha l[1] l[2])                   

                             (1/2) /    1
  (h[1] beta alpha l[1] l[2])      |d - -
                                   \    2

     1  (1/2)  /   / (1/2)   \      / (1/2)   \\ \\\     
   + - 3      { tan\3      xi/ + sec\3      xi/ }|||^3 + 
     2         \                               / //|     
                                                   /     

              /                                                 
       1      |                                           (1/2) 
  ----------- \alpha l[1] l[2] (h[1] beta alpha l[1] l[2])      
  2 beta h[1]                                                   

          /                 /                  2\
   (1/2) {      / (1/2)   \ |       / (1/2)   \ |
  3       \6 tan\3      xi/ \1 + tan\3      xi/ /

                                    2
          / (1/2)   \    / (1/2)   \ 
   + 3 sec\3      xi/ tan\3      xi/ 

                      /                  2\\ \
          / (1/2)   \ |       / (1/2)   \ | }|
   + 3 sec\3      xi/ \1 + tan\3      xi/ // /
simplify(%);
Error, (in simplify/power) invalid input: ln expects its 1st argument, x, to be of type algebraic, but received {(sin(3^(1/2)*xi)+1)/cos(3^(1/2)*xi)}
 

How to "stack" the procedures generated by dsolve ...

Asked by: mmcdara 7891
ode parameters dsolve
September 09 2019
0  3

Hi, 

I solve numerically an ode for different values of its parameters ( dsolve(..., numeric, parameters=[...] ) and I would like to "stack" the different solutions in a container (the container (list, vector, table) is of no matter).

Here is a notional example where I try to construct a sequence where the first element should be the solution when the parameter is equal to 1 and second one when the parameter is equal to 2.
It happens that some "premature evaluation" seems to occur which makes the two elements identical.

Please do not pay attention to the obvious simplicity of the problem: the true one is more complicated but can be illustrated by the on below.

Thanks in advance
 

> 

restart:
> 

f := dsolve({diff(x(t),t)=A*t, x(0)=0}, numeric, parameters=[A]);

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := [A = A]; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 24, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 1, (2) = 1, (3) = 0, (4) = 0, (5) = 1, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 0, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = .0, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..2, {(1) = 0., (2) = Float(undefined)})), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..1, {(1) = .1}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8]), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..1, {(1, 1) = .0, (2, 0) = .0, (2, 1) = .0, (3, 0) = .0, (3, 1) = .0, (4, 0) = .0, (4, 1) = .0, (5, 0) = .0, (5, 1) = .0, (6, 0) = .0, (6, 1) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t)]`; YP[1] := Y[2]*X; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t)]`; YP[1] := Y[2]*X; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0) ] )) ] ); _y0 := Array(0..2, {(1) = 0., (2) = 0.}); _vmap := array( 1 .. 1, [( 1 ) = (1) ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; _dat[4][26] := _EnvDSNumericSaveDigits; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, x(t)], (4) = [A = A]}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

 (1)
> 

k := 1:
for a in [1, 2] do
  f(parameters=[a]):
  printf(" f(1) = %a\n", f(1)):
  g||k := unapply(f(t), t):
  for kk from 1 to k do
    printf("g%d(1) = %a\n", kk, g||kk(1)):
  end do:
  k := k+1:
  print():
end do:

 f(1) = [t = 1., x(t) = .500000000000001]
g1(1) = [t = 1., x(t) = .500000000000001]

  

 

 f(1) = [t = 1., x(t) = .999999999999999]
g1(1) = [t = 1., x(t) = .999999999999999]
g2(1) = [t = 1., x(t) = .999999999999999]

  

(2)
> 

# how must I correct this in order to prevent the
# "over writting" of g1 when g2 is instanciated
# and get
#
#  f(1) = [t = 1., x(t) = .999999999999999]
# g1(1) = [t = 1., x(t) = 0.5]
# g1(2) = [t = 1., x(t) = .999999999999999]
#
> 

 


 

Download Stack_procedures.mw

Bug using Heaviside and PDF in statistics package...

Asked by: markweitzman 110
integration heaviside statistics
September 08 2019
1  2

Integrating a positive definite function (normal distribution) and a >= 0 function (Heaviside) should not return a negative value.

 

 

with(Statistics);
X := RandomVariable(Normal(1, sqrt(2.25)));
int(PDF(X, x)*Heaviside(x^7-5*x^4-3*x+1), x = -infinity .. infinity);
                         -0.08507120131

 

 

Bug.mw

How to Generate a Longest Path and Longest Cycle i...

Asked by: DJJerome1976 540
graph cycle path
September 08 2019
1  3

I'm doing some work with undirected and directed graphs, and I'm needing to find examples of longest paths and longest cycles.  How may this be done in Maple?

How to Solve a system of coupled ODEs in Maple ?...

Asked by: Mayo_Kun 5
ode dsolve numeric
September 08 2019
1  5

Hello,

I am trying to solve a set of coupled ODEs in the following code (ODE_Prob.mw). But I am getting an error.

Please, what does this error mean and what is the solution.

Thanks in advance. 
 

> 

restart; PDEtools[declare](f(eta), prime = eta, theta(eta), prime = eta); Digits := 16; inf := 9

f(eta)*`will now be displayed as`*f

  

`derivatives with respect to`*eta*`of functions of one variable will now be displayed with '`

  

theta(eta)*`will now be displayed as`*theta

 (1)
> 

equ1 := (2*n+1)*f(eta)*(diff(theta(eta), eta))/(2*(n+1)) = (diff(theta(eta), eta, eta))/Pr^(2/(n+1))

(2*n+1)*f(eta)*(diff(theta(eta), eta))/(2*n+2) = (diff(diff(theta(eta), eta), eta))/Pr^(2/(n+1))

 (2)
> 

equ2 := (1/2)*(diff(f(eta), eta))^2-(2*n+1)*f(eta)*(diff(f(eta), eta, eta))/(2*n+2) = -n*(diff(f(eta), eta, eta))^(n-1)*(diff(f(eta), eta, eta, eta))+theta(eta)

(1/2)*(diff(f(eta), eta))^2-(2*n+1)*f(eta)*(diff(diff(f(eta), eta), eta))/(2*n+2) = -n*(diff(diff(f(eta), eta), eta))^(n-1)*(diff(diff(diff(f(eta), eta), eta), eta))+theta(eta)

 (3)
> 

indets(equ1);

{Pr, eta, n, Pr^(2/(n+1)), diff(diff(theta(eta), eta), eta), diff(theta(eta), eta), f(eta), theta(eta)}

  

{eta, n, (diff(diff(f(eta), eta), eta))^(n-1), diff(diff(diff(f(eta), eta), eta), eta), diff(diff(f(eta), eta), eta), diff(f(eta), eta), f(eta), theta(eta)}

 (4)
> 

Bcs := f(0) = 0, (D(f))(0) = 0, (D(f))(inf) = 1, theta(0) = 1, theta(inf) = 0;

f(0) = 0, (D(f))(0) = 0, (D(f))(9) = 1, theta(0) = 1, theta(9) = 0

  

1.301763965

  

1.4

 (5)
> 

equ1;

.7916666666666667*f(eta)*(diff(theta(eta), eta)) = .8027059130737796*(diff(diff(theta(eta), eta), eta))

  

(1/2)*(diff(f(eta), eta))^2-.7916666666666667*f(eta)*(diff(diff(f(eta), eta), eta)) = -1.4*(diff(diff(f(eta), eta), eta))^.4*(diff(diff(diff(f(eta), eta), eta), eta))+theta(eta)

 (6)
> 

SolP1 := dsolve({Bcs, equ1, equ2}, numeric); SolP1(0); SolP1(inf)

Error, (in dsolve/numeric/BVPSolve) unable to store '-.4744942626692375+1.460343180701223*I' when datatype=sfloat

  

SolP1(0)

  

SolP1(9)

 (7)
> 

p1 := plots:-odeplot(SolP1, [eta, f(eta)], 0 .. inf); p2 := plots:-odeplot(SolP1, [eta, theta(eta)], 0 .. inf); p3 := plots:-odeplot(SolP1, [eta, diff(f(eta), eta)], 0 .. inf); plots:-display(Array([p1, p2, p3]))

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

  

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

  

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

  

Error, (in plots:-display) element 1 of the rtable is not a valid plot structure

  
> 

``


 

Download ODE_Prob.mw

 

First 643 644 645 646 647 648 649 Last Page 645 of 2425