1、二阶常微分方程的数值求解,一. 教学要求,掌握利用降阶把二阶常微分方程转化为一阶微分 方程组，再利用Euler方法数值求解,并能利用MATLAB 软件进行数值计算和符号运算。,二. 教学过程,考虑如下的二阶微分方程初值问题,利用Euler方法求解上述方程组可得如下数值格式,利用四阶R-K方法求解上述方程组可得如下数值格式,例1：用 Euler 法求解如下初值问题,当 h=0.1，即 n=20 时，Matlab 源程序见 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
2、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_sys1.m,数值解与真解如下图,例2：利用4阶R-K方法求解例1，并与Euler方法进行比较。,解 当 h=0.1，即 n=20 时，R-K方法的Matlab 源程序见 RK_sys1.m，数值结果见下图,function w=rightf_sys1(x,y,z)w=y;,clc;clear; h=0.1; a=0;b=2; x=a:h:b; Euler_y(1)=
3、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); %* R-K4 Method*%K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i); % K1 and L1K2=RK_z(i)+0.5*h*L1;,rightf_sys1.m,RK_sys1.m,L2=rightf_sys1(x(i
4、)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1); % K2 and L2K3=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 L3K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); % K4 and L4RK_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*
5、L2+2*L3+L4); end plot(x,Euler_y,r+,x,exp(-x),k-,x,RK_y,b*); xlabel(Variable x); ylabel(Variable y);,例3：分别用 Euler 法和R-K4求解如下初值问题,解：,当 h=0.1，即 n=20 时，Matlab 源程序见 RK_sys2.m, 数值结果如下图,function w=rightf_sys2(x,y,z) w=-y+2*exp(-x)*(x-1);,clc;clear; h=0.1; a=0;b=2; x=a:h:b; Euler_y(1)=1; Euler_z(1)=1; RK_y(
6、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,rightf_sys1.m,RK_sys2.m,K2=RK_z(i)+0.5*h*L1; L2=rightf_sys2(x(i)+
7、0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1); % K2 and L2K3=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 L3K4=RK_z(i)+h*L3; L4=rightf_sys2(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); % K4 and L4RK_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
8、+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);,dsolve 的调用格式,y=dsolve(eq1,eq2, . ,cond1,cond2, . ,v),其中 y 为输出， eq1、eq2、.为微分方程，cond1、cond2、.为初值条件，v 为自变量，如果不指定v作为自变量，则默认t为自变量。,例 4：求微分方程 的通解，并验证。, y=dsolve(Dy+2*x*y=x*exp(-x2),x), syms x; diff(
9、y)+2*x*y - x*exp(-x2),利用dsolve 函数求微分方程解析解,几点说明,如果省略初值条件，则表示求通解；,如果省略自变量，则默认自变量为 t,dsolve(Dy=2*x,x); dy/dx = 2x dsolve(Dy=2*x); dy/dt = 2x,若找不到解析解，则返回其积分形式。,微分方程中用 D 表示对 自变量 的导数，如：,Dy y； D2y y； D3y y,例 5：求微分方程 在初值条件 下的特解，并画出解函数的图形。, y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x) ezplot(y);,例 6,在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),注意：只有很少一部分微分方程（组）能求出解析解。 大部分微分方程（组）只能利用数值方法求数值解。,作业,利用Euler方法和R-K方法求解一个二阶常微分初值问题，并比较数值结果，计算数值解和解析解的误差。 利用dsolve函数求解一些微分方程的通解,
