二阶常微分方程数值求解

解：令z y',则该初值问题可以转化为
y'(x) z(x),
x[0,2]
z'(x) y(x),
x[0,2]
z(0) 1
z0,
y(0) 1.
当 h=0.1，即 n=20 时，Matlab 源程序见 Euler_sys1.m
Euler_sys1.m
clc;clear; h=0.1; a=0;b=2; x=a:h:b; y(1)=1; z(1)=-1; for i=1:length(x)-1
z(0) 1 z0, y(0) 1.
当 h=0.1，即 n=20 时，Matlab 源程序见 RK_sys2.m, 数值结
果如下图
rightf_sys1.m
function w=rightf_sys2(x,y,z) w=-y+2*exp(-x)*(x-1);
RK_sys2.m
clc;clear; h=0.1; a=0;b=2; x=a:h:b; Euler_y(1)=1; Euler_z(1)=1; RK_y(1)=1; RK_z(1)=1; for i=1:length(x)-1 %**** Euler Method ****%
zk
yk yk1 hzk , zk1 hf ( xk , yk , zk ),
k
1, 2
xk xk1 h.
其中yk是y(xk )的近似，zk是y '(xk )的近似
利用四阶R-K方法求解上述方程组可得如下 数值格式
h yk1 yk 6 ( K1 2K 2 2K 3 K 4 ),
K 4 zk hL3 , L4 f ( xk h, yk hK 3 , zk hL3 ).
其中yk是y(xk )的近似，zk是y '(xk )的近似
k 1, 2
例1：用 Euler 法求解如下初值问题
d2y
dx2 y
y(0)1, y'(0)1
该问题的真解为 y e x
x[0,2]
Euler_y(i+1)=Euler_y(i)+h*Euler_z(i); Euler_z(i+1)=Euler_z(i)+h*rightf_sys2(x(i),Euler_y(i),Euler_z(i)); %***** R-K4 Method*****% K1=RK_z(i); L1=rightf_sys2(x(i),RK_y(i),RK_z(i)); % K1 and L1
zk 1
zk
h 6
( L1
2 L2
2 L3
பைடு நூலகம்
L4
),
K1 zk , L1 f ( xk , yk , zk ),
K
2
zk
h 2
L1 ,
h
h
h
L2 f ( xk 2 , yk 2 K1 , zk 2 L1 ),
K3
zk
h 2
L2 ,
h
h
h
L3 f ( xk 2 , yk 2 K 2 , zk 2 L2 ),
% K1 and L1
L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1); % K2 and L2
K3=RK_z(i)+0.5*h*L2; L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2); % K3 and L3
y(i+1)=y(i)+h*z(i); z(i+1)=z(i)+h*y(i); end plot(x,y,'r+',x,exp(-x),'k-'); xlabel('Variable x'); ylabel('Variable y');
数值解与真解如下图
例2：利用4阶R-K方法求解例1，并与Euler方法 进行比较。
例3：分别用 Euler 法和R-K4求解如下初值问题
d dx 2y 2 y2ex(x1) y(0)1, y'(0)1
x[0,2]
该问题的真解为 y cos x xe x
解：令z y',则该初值问题可以转化为
y'(x) z(x),
x[0,2]
z'(x) y(x)2ex(x1), x[0,2]
二阶常微分方程的数值求解
一. 教学要求
掌握利用降阶把二阶常微分方程转化为一阶微分 方程组，再利用Euler方法数值求解,并能利用MATLAB 软件进行数值计算和符号运算。
二. 教学过程
考虑如下的二阶微分方程初值问题
d 2 y d x 2 f(x ,y ,y '),y (a ) y 0,y '(a ) y 1 ,x [a ,b ]
%**** Euler Method ****% Euler_y(i+1)=Euler_y(i)+h*Euler_z(i); Euler_z(i+1)=Euler_z(i)+h*Euler_y(i);
%***** R-K4 Method*****% K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); K2=RK_z(i)+0.5*h*L1;
K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); % K4 and L4
RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4); RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4); end plot(x,Euler_y,'r+',x,exp(-x),'k-',x,RK_y,'b*'); xlabel('Variable x'); ylabel('Variable y');
解 当 h=0.1，即 n=20 时，R-K方法的Matlab 源程序见
RK_sys1.m，数值结果见下图
rightf_sys1.m
function w=rightf_sys1(x,y,z) w=y;
RK_sys1.m
clc;clear; h=0.1; a=0;b=2; x=a:h:b; Euler_y(1)=1; Euler_z(1)=-1; %初值 RK_y(1)=1; RK_z(1)=-1; %初值 for i=1:length(x)-1
若令 z y ' , 则上述初值问题可以转化为如下一
阶微分方程组初值问题
y'(x)z(x),
x[a,b]
z'(x) f(x, y(x),z(x)), x[a,b]
z(a)y1 z0, y(a)y0
利用Euler方法求解上述方程组可得如下数 值格式
y(a) y0 , y '(a) z0