ode .

“msharpmath”, The Simple is the Best

[100] 006 Tutorial ode              ordinary differential equations

Click the above button for full-text tutorial of   ‘ode’  or   ‘ivp’.

//==================================================================================
//  Umbrella  ‘ode’  for Initial Value Problems (IVP)
//==================================================================================
DE   :  y^(n) = F(t,y,y’,y”, … , y^(n-1)),                               n <= 4  up to 4-th order ODE
ICs  :  y = yo,   y’ = yo_1,   y” = yo_2,   …   ,   y^(n-1) = yo_(n-1)

ode .t[n=101,g=1](a,b) ( << opt >>, DE, ICs, DE, ICs, … ) .Spoke.Spoke..

ode  .t[n=101,g=1](a,b) ( << opt >>, DE, DE, … ICs, ICs, … ) .Spoke.Spoke..

//————————————————————————————————————
//  Spokes
//————————————————————————————————————
.plot(functions of t,y,y’,…)
.return(functions of t,y,y’,…)       extract the solution and its derivatives
.togo ( X=x, … ) extract the specified field

.euler      //  .case(1)  1st-order Euler Method
.euler2     //  .case(2)  2nd-order modified Euler method
.pc2        //  .case(3)  2nd-order predictor-corrector
.rs2        //  .case(4)  2nd-order Ralston method
.rk3        //  .case(5)  3rd-order Runge-Kutta method
.r4         //  .case(6)  one of 4th-order Runge-Kutta methods
.k4         //  .case(7)  one of 4th-order Runge-Kutta methods
.rkg4     //  .case(8)  4th-order Runge-Kutta-Gill method
.rkf4/rk //  .case(9)  4th-order Runge-Kutta-Fehlberg method
.rkf5       //  .case(10) 5th-order Runge-Kutta-Fehlberg method (no adaptation)

.ab(n)      //  .case(11~16)  nth-order Adams-Bashforth method
.am(n)      //  .case(21~26)  nth-order Adams-Molton method

.apc        //  .case31)  Adams-predictor-corrector method
.mpc        //  .case(32)  Milne-predictor-corrector method
.hpc        //  .case(33)  Hamming-predictor-corrector method
.mopup      //  .case(34)  Adams-Molton with mop-up method

.stopif( condition)    // terminate current step if condition is true
.update(expr,expr,…)   // update after each step of integration

.stiff(tol=0.01)     // applied to stiff ODEs, variable step size

//————————————————————————————————————
//  Implicit  1st-order ODE : Syntax
//————————————————————————————————————

ode# .t[n=101,g=1](a,b) ( DE, ICs ) .Spoke.Spoke..

//————————————————————————————————————
//  Examples
//————————————————————————————————————

%> Euler and modified Euler ===================================
#> ode.x[ 3](0,0.2) ( y’ = y,  y = 1 ).euler2.return(x,y,-y+exp(x)); // method(2)
#> ode.x[11](0,0.2) ( y’ = y,  y = 1 ).euler .return(x,y,-y+exp(x)); // method(1)
ans =
[             0             1             0 ]
[           0.1         1.105    0.00017092 ]
[           0.2         1.221    0.00037776 ]
ans =
[             0             1             0 ]
[          0.02          1.02    0.00020134 ]
[          0.04        1.0404    0.00041077 ]
[          0.06        1.0612    0.00062855 ]
[          0.08        1.0824    0.00085491 ]
[           0.1        1.1041     0.0010901 ]
[          0.12        1.1262     0.0013344 ]
[          0.14        1.1487     0.0015881 ]
[          0.16        1.1717     0.0018515 ]
[          0.18        1.1951     0.0021248 ]
[           0.2         1.219     0.0024083 ]
Even with one-tenth of step size, the modified-Euler method is more accurate than the Euler method, as expected.

%> severe variation ===================================
#> ode.x[51](0,5) ( y’ = exp(x)*y,  y = exp(1) )
.return( x,y, exp(exp(x)), y/exp(exp(x)) );
ans =
[             0        2.7183        2.7183             1 ]
[           0.1        3.0197        3.0197             1 ]
[           0.2        3.3919        3.3919             1 ]

[           4.8   5.9058e+052   5.9061e+052       0.99994 ]
[           4.9   2.0955e+058   2.0956e+058       0.99993 ]
[             5   2.8509e+064   2.8511e+064       0.99991 ]
Note 2.85e+64 at the end where fnum/fex is still near 1.0.

%> Van der Pol equation ===================================
#> mu=1;
#> ode .t(0,20) ( y” = mu*(1-y*y)*y’ – y, y=2, y’=0 )
.plot( y,y’ )
.return(t,y,y’);

    Van der Pol (non-stiff only)

%> plot and data-extraction ===================================
#{
ode .x[11](0,pi) ( y” = -y,  y = 0, y’ = 1 )
.plot( y,y’,-y+sin(x) )
.return( x,y,y’, y*y+y’*y’-1 );
#}
ans =
[             0             0             1             0 ]
[      0.314159      0.309021      0.951058   4.9971e-006 ]
[      0.628319      0.587794      0.809017  9.99423e-006 ]
[      0.942478      0.809029      0.587781  1.49914e-005 ]
[       1.25664       0.95107      0.309007  1.99886e-005 ]
[        1.5708       1.00001 -1.68818e-005  2.49858e-005 ]
[       1.88496      0.951065     -0.309041   2.9983e-005 ]
[       2.19911      0.809017     -0.587815  3.49803e-005 ]
[       2.51327      0.587775     -0.809049  3.99775e-005 ]
[       2.82743      0.308995     -0.951087  4.49748e-005 ]
[       3.14159 -3.37641e-005      -1.00002  4.99722e-005 ]

   plot of  y,y’, y*y+y’*y’-1

%> Coupled ODE =========================================


#> ode.x[21](0,1)
(
u” = -v’-18*u + 56/3*cos(x)+18*exp(x),  u=3, u’=4,
v” = -4*v-5*u’ – 10*sin(2*x),           v=2, v’=-40/3
)
.plot( u,v, u’,v’,
-u + ( cos(x)+cos(2*x)+sin(3*x)+exp(x) ),
-v + ( 5/3*sin(x)-7*sin(2*x)+3*cos(3*x)-exp(x) ) )
.return( u,u’,v,v’ );

 

%> Coupled ODE =========================================
//  3x” – 2y” = -2x’ + y’ + x – 6y – t^2 + 10t + 5
//  -4x” + 3y” = 3x’ – 5y’ – x + 2y + t^2 – 8t – 3,
//  x = 1, x’ = -1, // x = exp(-t) + t^2
//  y = 1, y’ = -1  // y = exp(-2t) + t

The above coupled ODE is solved by elimination as follows.

eps = 0.002;;
A = [ 3,-2; -4,3 ];

ode .t(0,1) (
<<
f = -2*x’ + y’ + x – 6*y – t^2 + 10*t + 5, // for 2nd-order
g =  3*x’ – 5*y’ – x + 2*y + t^2 – 8*t – 3, // x’, y’ are here
B = A \ [ f ; g ]
>>,
x” = B(1),
y” = B(2), x = 1, y = 1, x’ = -1, y’ = -1
)
.plot( x,y,
exp(-t) + t^2 + eps,
exp(-2*t) + t + eps );

%> phase diagram =========================================
#> .hold;
ode  .x(0,10)        ( y’ = y*exp(-0.1*x)*sin(x), y = 1   ).plot(y);
ode  .x(0,10)        ( y’ = y*exp(-0.1*x)*sin(x), y = 0.5 ).plot+(y);
plot .x[21](0,10).y[21](0,6) ( 1, y*exp(-0.1*x)*sin(x)  ).phase[0.2,3];
#> plot;

%> phase diagram of Lorentz equation ============================
#{
ode.t[5001](0,50) (
x’ = 10*(y-x),
y’ = 28*x-y-x*z,
z’ = x*y-8/3*z,

x = 0, y = 1, z = 0
).plot@(x,z);
#}

%> phase diagram ==========================================
#{
t = (0,100).span(1001);;
(a,b,c) = (0.2, 0.2, 5);;
ode[t] (
x’ = -y-z,      x = 1,
y’ = x+a*y,     y = 1,
z’ = b+z*(x-c), z = 1
).plot@(x+1,y+1,z+1);
#}
#{
(a,b,c) = (0.2, 0.2, 2.5);;
ode[t] (
x’ = -y-z,      x = 1,
y’ = x+a*y,     y = 1,
z’ = b+z*(x-c), z = 1
).plot+@(x+20,y,z);
#}

%> bouncing motion of a ball ===================================
#> to = 0;;
#> (yo,vo) = (0, 20);;
#> ntouch = 0;;
.hold;
#> plot.circle(0.1).thick(9);; // ball

#{
while( ++ntouch <= 10 ) {
ode .t(to,to+10) ( y” = -9.8, y = yo, y’ = vo )
//         .stopif( (y’ < 0 && y <= 0) )
.stopif( (y’ < 0 && y <= 0)  || (y’ > 0 && y >= 10) )
.update( to = t, yo = y, vo = -0.9*y’ )
.plot+(y);
}
#}
plot;

%> fox and rabbit ===========================================
#{
double u(t) = sqrt(1+t)*cos(t);
double v(t) = sqrt(1+t)*sin(t);
eps = 1.e-6;
del = 100;  // assume initial distance far away
ntau = 1001;  // accurate-sensitive, thus many steps
#}

#{
.hold;
ode.t[ntau](0,6) (
<<
rx = u(t), // rabbit’s x-coordinate
ry = v(t), // rabbit’s y-coordinate
d = .max(eps, .norm2(rx-x,ry-y)), // distance
s = (1.1/2)* .norm2( rx/(1+t)-2*ry, ry/(1+t)+2*rx )
>>,
x’ = s/d * (rx-x),
y’ = s/d * (ry-y),
x = 3, y = 0  // fox’s initial position
)
.stopif( d > del ) // d’ > 0 condition is realized
.update( del = d, tau = t, xtau = x, ytau = y ) // make a new distance
.plot@(x,y);
#}

#{
plot .@t(0,tau) ( u(t),v(t) );
plot .circle(0.1).xmove(3); // fox
plot .circle(0.1).xmove(1); // rabbit
plot;

tau; del; xtau; ytau;
#}

tau =         5.06475
del =    0.0014870192
xtau =       0.8485762
ytau =      -2.3118532

%> implicit ODE ===========================================


A = ode# .t(1,10) (
t*y*y*(y’)^3 – y^3*(y’)^2+t*(t*t+1)*(y’)-t*t*y = 0,
y = sqrt(3/2)
)
.guess(1)
.plot( y, sqrt(t*t+0.5) )
.return( t,y, -y+sqrt(t*t+0.5) );

#> A ;

A =
[             1       1.22474             0 ]
[          1.09       1.30037   -0.00109735 ]
[          1.18       1.37777   -0.00212642 ]

[          9.82       9.90325    -0.0578221 ]
[          9.91       9.99355    -0.0583569 ]
[            10       10.0839    -0.0588916  ]

===========================================

Stiff ODEs

examples are for the van der Pol equations

 

#> ode.x(0,20) ( y” = 1*(1-y^2)*y’ – y, y = 2, y’ = 0 ).plot;

Figure 11a  Van der Pol equation with   mu=1

#> ode.x[1000](0,60) ( y” = 10*(1-y^2)*y’ – y, y = 2, y’ = 0 ).stiff(0.01) .return(x,y).skip(10).plot;

Figure 11b  Van der Pol equation with   mu=10

#> ode.x[15000](0,300) ( y” = 100*(1-y^2)*y’ – y, y = 2, y’ = 0 ).stiff .return(x,y).skip(10).plot;

Figure 11c  Van der Pol equation with   mu=100

Comments are closed.