二阶常微分方程的数值求解

dy x2 2 xy xe 例 4：求微分方程 的通解，并验证。 dx
>> syms x; diff(y)+2*x*y - x*exp(-x^2) >> y=dsolve('Dy+2*x*y=x*exp(-x^2)','x')
几点说明
微分方程中用 D 表示对 自变量 的导数，如： Dy y'； D2y y''； D3y y'''
当 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 ****% 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
例3：分别用 Euler 法和R-K4求解如下初值问题
d2 y 2 y 2e x ( x 1) dx y(0) 1, y '(0) 1
该问题的真解为 y cos x xe x
x [0, 2]
解：
令z y ', 则该初值问题可以转化为 y '( x ) z ( x ), x [0, 2] z '( x ) y( x ) 2e x ( x 1), x [0, 2] z (0) 1 z0 , y(0) 1.
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 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');
利用Euler方法求解上述方程组可得如下数 值格式
zk y(a ) y0 , y '(a ) z0 yk yk 1 hzk , zk 1 hf ( xk , yk , zk ), k 1, 2 xk xk 1 h.
其中yk 是y( xk )的近似，zk 是y '( xk )的近似
利用四阶R-K方法求解上述方程组可得如下 数值格式
h yk 1 yk ( K1 2 K 2 2 K 3 K 4 ), 6 h zk 1 zk ( L1 2 L2 2 L3 L4 ), 6 K1 zk , L1 f ( xk , yk , zk ), K 2 zk h L1 , L2 f ( xk h , yk h K1 , zk h L1 ), 2 2 2 2 h h h h K 3 zk L2 , L3 f ( xk , yk K 2 , zk L2 ), 2 2 2 2 K 4 zk hL3 , L4 f ( xk h, yk hK 3 , z k hL3 ).
K2=RK_z(i)+0.5*h*L1; L2=rightf_sys2(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_sys2(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2); % K3 and L3 K4=RK_z(i)+h*L3; L4=rightf_sys2(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,cos(x)+x.*exp(-x),'k-',x,RK_y,'b*'); xlabel('Variable x'); ylabel('Variable y');
%***** 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;
% 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 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');
二阶常微分方程的数值求解
一. 教学要求
掌握利用降阶把二阶常微分方程转化为一阶微分 方程组，再利用Euler方法数值求解,并能利用MATLAB 软件进行数值计算和符号运算。
二. 教学过程
考虑如下的二阶微分方程初值问题
d2 y f ( x , y , y ') , y(a ) y0 , y '(a ) y1 , x [a , b] 2 dx
数值解与真解如下图
例2：利用4阶R-K方法求解例1，并与Euler方法 进行比较。
解 当 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 %**** Euler Method ****% Euler_y(i+1)=Euler_y(i)+h*Euler_z(i); Euler_z(i+1)=Euler_z(i)+h*Euler_y(i);
若令
z y'
, 则上述初值问题可以转化为如下一
阶微分方程组初值问题
x [a , b] y '( x ) z ( x ), z '( x ) f ( x , y( x ), z( x )), x [a , b] z (a ) y z , y(a ) y 1 0 0
如果省略初值条件，则表示求通解； 如果省略自变量，则默认自变量为
dsolve('Dy=2*x','x'); dsolve('Dy=2*x');
t
％ dy/dx = 2x ％ dy/dt = 2x
若找不到解析解，则返回其积分形式。
例 5：求微分方程 xy ' y e x 0 在初值条件 y(1) 2e
下的特解，并画出解函数的图形。
>> y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') >> ezplot(y);
例6
d2 y 2 cos( 2 x ) y 求二阶常微分方程 dx 的通解 y(0) 1, y '(0) 0
作业 利用Euler方法和R-K方法求解一个 二阶常微分初值问题，并比较数 值结果，计算数值解和解析解的 误差。 利用dsolve函数求解一些微分方程 的通解
k 1, 2
其中yk 是y( xk )的近似，zk 是y '( xk )的近似
例1：用 Euler 法求解如下初值问题
d2 y y 2 dx y(0) 1, y '(0) 1
该问题的真解为 y e x
x [0, 2]
解： 令z y ', 则该初值问题可以转化为
x [0, 2] y '( x ) z ( x ), x [0, 2] z '( x ) y( x ), z (0) 1 z , y(0) 1. 0
当 h=0.1，即 n=20 时，Matlab 源程序见 Euler_sys1.m
Euler_sys1.m
源自文库
在Matlab中的命令窗口中输入下面的命令
>> syms x y >> S=dsolve('D2y=cos(2*x)-y','y(0)=1','Dy(0)=0','x')
则可以得到如下的结果
S=
4/3*cos(x)-1/3*cos(2*x)
注意：只有很少一部分微分方程（组）能求出解析解。 大部分微分方程（组）只能利用数值方法求数值解。
利用dsolve 函数求微分方程解析解
dsolve 的调用格式
y=dsolve('eq1','eq2', ... ,'cond1','cond2', ... ,'v') 其中 y 为输出， eq1、eq2、...为微分方程，cond1、 cond2、...为初值条件，v 为自变量，如果不指定v作为自变 量，则默认t为自变量。