수치해석68🧮 7-6 msharpmath 활용 (상미분방정식)

[@ptostudy]

동영상에 포함된 명령어 모음.
특별히 이 댓글에는 상미분방정식을 푸는 서브루틴을 포함하고 있어서
상당히 길지만 독자가 수정해볼 수 있습니다.
복사 후 편집창에 붙여넣기하고 Fn key 누른 상태에서 F5 key 누르면 실행됩니다.
#> ode .t[1+100](0,3) ( y' = -y-5*exp(-t)*sin(5*t), y = 1 )
.plot( y, y/(exp(-t)*cos(5*t)) )
.return( y-exp(-t)*cos(5*t) ).maxentry;
// prototype
matrix user_ode (double t,tf,h, matrix y,fun(double t,matrix y));
matrix adaptive_ode (double t,tf,h, matrix y,fun(double t,matrix y), double tol);
//========================================
// user_odecoef::
//========================================
double user_odecoef(double icase, matrix &A,&B,&C) {
switch( icase ) {
case 1: // Euler = 1st RK
A = [1]; B = [0];
C = [ 0 ];
break;
case 2: // (modified Euler) predictor-corrector = 2nd RK
A = [1,1]/2; B = [0,1];
C = [ 0; 1,0 ];
break;
case 3: // polygon = 2nd RK
A = [0,1]; B = [0,1]/2;
C = [ 0; 0.5,0 ];
break;
case 4: // Ralston = 2nd RK
A = [1,2]/3; B = [0,3]/4;
C = [ 0; 0.75,0 ];
break;
case 5: // 3rd RK
A = [1,4,1]/6; B = [0,1,2]/2;
C = [ 0; 0.5; -1,2,0 ];
break;
case 6: // Runge = 4th RK
A = [1,2,2,1]/6; B = [0,1,1,2]/2;
C = [ 0; 0.5; 0,0.5; 0,0,1,0 ];
break;
case 7: // Kutta = 4th RK
A = [1,3,3,1]/8; B = [0,1,2,3]/3;
C = [ 0; 1/3; -1/3,1; 1,-1,1,0 ];
break;
case 8: // Gill = 4th RK
s = 1/sqrt(2);
A = [0.5,1-s,1+s,0.5]/3; B = [0,1,1,2]/2;
C = [ 0; 0.5; -0.5+s, 1-s; 0, -s, 1+s, 0 ];
break;
case 9: // Fehlberg = 4th RK
A = [ 25/216, 0, 1408/2565, 2197/4104, -0.2 ];
B = [ 0, 0.25, 3/8, 12/13, 1 ];
C = [ 0;
0.25;
3/32, 9/32; // 3rd
1932/2197, -7200/2197, 7296/2197; // 4th
439/216, -8, 3680/513, -845/4104, 0 ]; // 5th
break;
case 10: // Fehlberg = 5th RK
A = [ 16/135, 0, 6656/12825, 28561/56430, -9/50, 2/55 ];
B = [ 0, 0.25, 3/8, 12/13, 1, 0.5 ];
C = [ 0;
0.25;
3/32, 9/32; // 3rd
1932/2197, -7200/2197, 7296/2197; // 4th
439/216, -8, 3680/513, -845/4104; // 5th
-8/27, 2, -3544/2565, 1859/4104, -11/40,0 ];
} 0;
return A.len; // A.len = A.length = A.longer
} 0;;
//========================================
// user_odeRK::
//========================================
// to : previous time at whicvh solutions are known
// h : magnitude of time step, h = t-to
// yo : solution at previous time step
// y : dependent variables for coupled ODEs, y(1),y(2), ... , y(numeq)
// A,B,C : coefficients for integration scheme
// K : generalized slopes for integration scheme
void user_odeRK( double to,h, // & means 'refer' not 'copy'
matrix &y,&yo, &A,&B,&C,&K, matrix fun(double t, matrix y)) {
K.row(1) = fun(to,yo); // B(1) = C.row(1) = 0
for.i(2,A.len) {
t = to + B(i)*h;
y = yo + h*( C.row(i) * K );
K.row(i) = fun(t,y);
} 0;
y = yo + h*( A * K );
} 0;;
//========================================
// user_ode:: for positive h
//========================================
// t : initial time
// tf : final time
// h : size of time step
// y : initial conditions, [ y(1), y(2), ... , y(n) ]
// fun : ODE in matrix form, [ y'(1)=f1, ... , y'(n)=fn ]
matrix user_ode(double t,tf,h, matrix y,fun(double t,matrix y)) {
icase = 9; // Fehlberg 4th, change this manually
sdim = user_odecoef(icase, A,B,C);
K = .zeros(sdim,y.len);
yo = fun(t,y);
if( yo.len != y.len ) "ode function is not valid in dim";
num = 0;
ysol = [ t, y ];
while(1) {
yo = y;
if( t+h >= tf ) h = tf-t ; // stay exact at the final time
user_odeRK(t,h, y,yo, A,B,C,K, fun);
t += h;
ysol _= [ t, y ]; // stack newly-obtained solutions

if( t >= tf || t.-.tf < 1.e-6 || ++num > 10000 ) break;
} 0;
return ysol;
} 0;;
//========================================
// adaptive_ode:: for positive h
//========================================
// t : initial time
// tf : final time
// h : size of time step
// y : initial conditions, [ y(1), y(2), ... , y(n) ]
// fun : ODE in matrix form, [ y'(1)=f1, ... , y'(n)=fn ]
matrix adaptive_ode(double t,tf,h, matrix y, fun(double t,matrix y), double tol) {
s = user_odecoef( 9, A5,B5,C5); K5 = .zeros(s,y.len);//4th
s = user_odecoef(10, A6,B6,C6); K6 = .zeros(s,y.len);
E6 = [ 1/360, -128/4275, -2197/75240, 0, 1/50, 2/55 ];
yo = fun(t,y);
if( yo.len != y.len ) "ode function is not valid in dim";
num = 0;
ysol = [ t, y, 0 ];
while(1) {
yo = y;
if( t+h > tf ) h = tf-t ; // stay exact at the final time
do {
user_odeRK(t,h, y,yo, A6,B6,C6,K6, fun);
yerr = h * |(E6*K6).maxentry|; // error estimate
if( yerr > tol) h *= 0.1; // reduce h where stiff
} while( yerr > tol );
// y = yo + h*(A5 * K6..(1:5)()); // if prefered 4th
t += h;
ysol _= [ t, y, yerr ]; // stack newly-obtained sol
if( yerr < tol ) h *= 1.25; // increase h where plain
if( t >= tf || t.-.tf < 1.e-6 || ++num > 100000 ) break;
if( .mod(num,1000) == 0 ) "adaptive_ode running" ;;
} 0;
return ysol;
} 0;;
// example 4
#> matrix.format("%13.8f");
#> double fex(x) = x+2*sqrt(3)-3* sqrt(2)*tanh(x/sqrt(2)+atanh(sqrt(2/3)));
#> x = (0,5) .step(0.1);;
#> ode.x[51](0,5) ( y''' = -1+y'^2, y = 0, y' = 0, y'' = 2/sqrt(3) )
.plot( y', y'' ) // note that y' and y''
.return( x, fex(x), y/fex(x) ) // changed to a matrix
.rowskip(-10) ; // ..(11:10:51)();
#> matrix fun4(double x, matrix y) = [ y(2), y(3), y(2)^2-1 ];
#> ysol = user_ode(0, 5, 0.1, [ 0, 0, 2/sqrt(3)], fun4 );
#> ysol.plot;
// example 5
#> matrix.format("%16.7e");
#> double ue(x) = cos(x)+cos(2*x)+sin(3*x)+exp(x);
#> double ve(x) = 5/3*sin(x)-7*sin(2*x)+3*cos(3*x)-exp(x);
#> x = (0,1) .step(0.1);;
#> ode.x[11](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 )
.return( x, u, u/ue(x), v, v/ve(x) ).rowskip(-2) ; // matrix
#> matrix fun5(double x, matrix y) =
[ y(3),y(4),-y(4)-18*y(1)+56/3*cos(x)+18*exp(x),
-4*y(2)-5*y(3)-10*sin(2*x) ];
#> ysol = user_ode(0,1, 0.1, [ 3,2,4,-40/3 ], fun5);;
#> ysol.plot;
// example 6
#> b = [ 30*cosd(45); 0; 30*sind(45) ];;
#> ode.t[5001](0,2) (
>,
theta' = B(1), theta = 0,
phi' = B(2), phi = 0,
psi' = B(3), psi = 0
) .plot(theta,phi,psi);
matrix fun6(double x, matrix y) {
A = [ 0, 1, sin(y(1)) ;
cos(y(2)), 0, -sin(y(2))*cos(y(1));
sin(y(2)), 0, cos(y(2))*cos(y(1)) ];
B = A \ [ 30*cosd(45); 0; 30*sind(45) ];
return B.tr;
} 0;;
#> ysol = user_ode( 0,2, 0.002, [ 0,0,0 ], fun6 );;
#> ysol. plot;
// example 7
#> 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).return(t,x,y,z).endrow;
matrix fun7(double x, matrix y) = [ 10*(y(2)-y(1)), 28*y(1)-y(2)-y(1)*y(3), y(1)*y(2)-8/3*y(3) ];
//matrix fun7(double x, matrix y) { return [ 10*(y(2)-y(1)), 28*y(1)-y(2)-y(1)*y(3), y(1)*y(2)-8/3*y(3) ]; } 0;;
#> ysol = user_ode( 0,50,0.01, [ 0,1,0 ], fun7 );;
#> ysol.plot; // need more work for parametrized plot
#> plot(ysol(*,2),ysol(*,3)); // A(*,j)=A.col(j)
#> plot(ysol.col(2),ysol.col(4));

[@ptostudy]

#> ode .t[1+100](0,3) ( y' = -y-5*exp(-t)*sin(5*t), y = 1 )
.plot( y, y/(exp(-t)*cos(5*t)) )
.return( y-exp(-t)*cos(5*t) ).maxentry;
// prototype
matrix user_ode (double t,tf,h, matrix y,fun(double t,matrix y));
matrix adaptive_ode (double t,tf,h, matrix y,fun(double t,matrix y), double tol);
//========================================
// user_odecoef::
//========================================
double user_odecoef(double icase, matrix &A,&B,&C) {
switch( icase ) {
case 1: // Euler = 1st RK
A = [1]; B = [0];
C = [ 0 ];
break;
case 2: // (modified Euler) predictor-corrector = 2nd RK
A = [1,1]/2; B = [0,1];
C = [ 0; 1,0 ];
break;
case 3: // polygon = 2nd RK
A = [0,1]; B = [0,1]/2;
C = [ 0; 0.5,0 ];
break;
case 4: // Ralston = 2nd RK
A = [1,2]/3; B = [0,3]/4;
C = [ 0; 0.75,0 ];
break;
case 5: // 3rd RK
A = [1,4,1]/6; B = [0,1,2]/2;
C = [ 0; 0.5; -1,2,0 ];
break;
case 6: // Runge = 4th RK
A = [1,2,2,1]/6; B = [0,1,1,2]/2;
C = [ 0; 0.5; 0,0.5; 0,0,1,0 ];
break;
case 7: // Kutta = 4th RK
A = [1,3,3,1]/8; B = [0,1,2,3]/3;
C = [ 0; 1/3; -1/3,1; 1,-1,1,0 ];
break;
case 8: // Gill = 4th RK
s = 1/sqrt(2);
A = [0.5,1-s,1+s,0.5]/3; B = [0,1,1,2]/2;
C = [ 0; 0.5; -0.5+s, 1-s; 0, -s, 1+s, 0 ];
break;
case 9: // Fehlberg = 4th RK
A = [ 25/216, 0, 1408/2565, 2197/4104, -0.2 ];
B = [ 0, 0.25, 3/8, 12/13, 1 ];
C = [ 0;
0.25;
3/32, 9/32; // 3rd
1932/2197, -7200/2197, 7296/2197; // 4th
439/216, -8, 3680/513, -845/4104, 0 ]; // 5th
break;
case 10: // Fehlberg = 5th RK
A = [ 16/135, 0, 6656/12825, 28561/56430, -9/50, 2/55 ];
B = [ 0, 0.25, 3/8, 12/13, 1, 0.5 ];
C = [ 0;
0.25;
3/32, 9/32; // 3rd
1932/2197, -7200/2197, 7296/2197; // 4th
439/216, -8, 3680/513, -845/4104; // 5th
-8/27, 2, -3544/2565, 1859/4104, -11/40,0 ];
} 0;
return A.len; // A.len = A.length = A.longer
} 0;;
//========================================
// user_odeRK::
//========================================
// to : previous time at whicvh solutions are known
// h : magnitude of time step, h = t-to
// yo : solution at previous time step
// y : dependent variables for coupled ODEs, y(1),y(2), ... , y(numeq)
// A,B,C : coefficients for integration scheme
// K : generalized slopes for integration scheme
void user_odeRK( double to,h, // & means 'refer' not 'copy'
matrix &y,&yo, &A,&B,&C,&K, matrix fun(double t, matrix y)) {
K.row(1) = fun(to,yo); // B(1) = C.row(1) = 0
for.i(2,A.len) {
t = to + B(i)*h;
y = yo + h*( C.row(i) * K );
K.row(i) = fun(t,y);
} 0;
y = yo + h*( A * K );
} 0;;
//========================================
// user_ode:: for positive h
//========================================
// t : initial time
// tf : final time
// h : size of time step
// y : initial conditions, [ y(1), y(2), ... , y(n) ]
// fun : ODE in matrix form, [ y'(1)=f1, ... , y'(n)=fn ]
matrix user_ode(double t,tf,h, matrix y,fun(double t,matrix y)) {
icase = 9; // Fehlberg 4th, change this manually
sdim = user_odecoef(icase, A,B,C);
K = .zeros(sdim,y.len);
yo = fun(t,y);
if( yo.len != y.len ) "ode function is not valid in dim";
num = 0;
ysol = [ t, y ];
while(1) {
yo = y;
if( t+h >= tf ) h = tf-t ; // stay exact at the final time
user_odeRK(t,h, y,yo, A,B,C,K, fun);
t += h;
ysol _= [ t, y ]; // stack newly-obtained solutions

if( t >= tf || t.-.tf < 1.e-6 || ++num > 10000 ) break;
} 0;
return ysol;
} 0;;
//========================================
// adaptive_ode:: for positive h
//========================================
// t : initial time
// tf : final time
// h : size of time step
// y : initial conditions, [ y(1), y(2), ... , y(n) ]
// fun : ODE in matrix form, [ y'(1)=f1, ... , y'(n)=fn ]
matrix adaptive_ode(double t,tf,h, matrix y, fun(double t,matrix y), double tol) {
s = user_odecoef( 9, A5,B5,C5); K5 = .zeros(s,y.len);//4th
s = user_odecoef(10, A6,B6,C6); K6 = .zeros(s,y.len);
E6 = [ 1/360, -128/4275, -2197/75240, 0, 1/50, 2/55 ];
yo = fun(t,y);
if( yo.len != y.len ) "ode function is not valid in dim";
num = 0;
ysol = [ t, y, 0 ];
while(1) {
yo = y;
if( t+h > tf ) h = tf-t ; // stay exact at the final time
do {
user_odeRK(t,h, y,yo, A6,B6,C6,K6, fun);
yerr = h * |(E6*K6).maxentry|; // error estimate
if( yerr > tol) h *= 0.1; // reduce h where stiff
} while( yerr > tol );
// y = yo + h*(A5 * K6..(1:5)()); // if prefered 4th
t += h;
ysol _= [ t, y, yerr ]; // stack newly-obtained sol
if( yerr < tol ) h *= 1.25; // increase h where plain
if( t >= tf || t.-.tf < 1.e-6 || ++num > 100000 ) break;
if( .mod(num,1000) == 0 ) "adaptive_ode running" ;;
} 0;
return ysol;
} 0;;
// example 4
#> matrix.format("%13.8f");
#> double fex(x) = x+2*sqrt(3)-3* sqrt(2)*tanh(x/sqrt(2)+atanh(sqrt(2/3)));
#> x = (0,5) .step(0.1);;
#> ode.x[51](0,5) ( y''' = -1+y'^2, y = 0, y' = 0, y'' = 2/sqrt(3) )
.plot( y', y'' ) // note that y' and y''
.return( x, fex(x), y/fex(x) ) // changed to a matrix
.rowskip(-10) ; // ..(11:10:51)();
#> matrix fun4(double x, matrix y) = [ y(2), y(3), y(2)^2-1 ];
#> ysol = user_ode(0, 5, 0.1, [ 0, 0, 2/sqrt(3)], fun4 );
#> ysol.plot;
// example 5
#> matrix.format("%16.7e");
#> double ue(x) = cos(x)+cos(2*x)+sin(3*x)+exp(x);
#> double ve(x) = 5/3*sin(x)-7*sin(2*x)+3*cos(3*x)-exp(x);
#> x = (0,1) .step(0.1);;
#> ode.x[11](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 )
.return( x, u, u/ue(x), v, v/ve(x) ).rowskip(-2) ; // matrix
#> matrix fun5(double x, matrix y) =
[ y(3),y(4),-y(4)-18*y(1)+56/3*cos(x)+18*exp(x),
-4*y(2)-5*y(3)-10*sin(2*x) ];
#> ysol = user_ode(0,1, 0.1, [ 3,2,4,-40/3 ], fun5);;
#> ysol.plot;
// example 6
#> b = [ 30*cosd(45); 0; 30*sind(45) ];;
#> ode.t[5001](0,2) (
>,
theta' = B(1), theta = 0,
phi' = B(2), phi = 0,
psi' = B(3), psi = 0
) .plot(theta,phi,psi);
matrix fun6(double x, matrix y) {
A = [ 0, 1, sin(y(1)) ;
cos(y(2)), 0, -sin(y(2))*cos(y(1));
sin(y(2)), 0, cos(y(2))*cos(y(1)) ];
B = A \ [ 30*cosd(45); 0; 30*sind(45) ];
return B.tr;
} 0;;
#> ysol = user_ode( 0,2, 0.002, [ 0,0,0 ], fun6 );;
#> ysol. plot;
// example 7
#> 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).return(t,x,y,z).endrow;
matrix fun7(double x, matrix y) = [ 10*(y(2)-y(1)), 28*y(1)-y(2)-y(1)*y(3), y(1)*y(2)-8/3*y(3) ];
//matrix fun7(double x, matrix y) { return [ 10*(y(2)-y(1)), 28*y(1)-y(2)-y(1)*y(3), y(1)*y(2)-8/3*y(3) ]; } 0;;
#> ysol = user_ode( 0,50,0.01, [ 0,1,0 ], fun7 );;
#> ysol.plot; // need more work for parametrized plot
#> plot(ysol(*,2),ysol(*,3)); // A(*,j)=A.col(j)
#> plot(ysol.col(2),ysol.col(4));