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

二二阶常微分方程的数常微分方程的数值求解求解一一.教学要求教学要求 掌握利用降阶把二阶常微分方程转化为一阶微分方程组，再利用Euler方法数值求解,并能利用MATLAB软件进行数值计算和符号运算。二二.教学过程教学过程q 考虑如下的二阶微分方程考虑如下的二阶微分方程初值问题初值问题q利用利用Euler方法求解上述方程组可得如下数方法求解上述方程组可得如下数值格式值格式q利用四阶利用四阶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 y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*y(i);endplot(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)=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 L1 K2=RK_z(i)+0.5*h*L1;rightf_sys1.mRK_sys1.mL2=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*K3RK_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);endplot(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(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=rightfsys2(x(i),RK_y(i),RK_z(i);%K1 and L1 rightf_sys1.mRK_sys2.mK2=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);endplot(x,Euler_y,r+,x,cos(x)+x.*exp(-x),k-,x,RK_y,b*);xlabel(Variable x);ylabel(Variable y);q 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(y)+2*x*y-x*exp(-x2)q利用利用dsolve dsolve 函数求微分方程解析解函数求微分方程解析解q 几点说明几点说明l 如果省略初值条件，则表示求通解；如果省略初值条件，则表示求通解；l 如果省略自变量，则默认自变量为如果省略自变量，则默认自变量为 t dsolve(Dy=2*x,x);dy/dx=2xdsolve(Dy=2*x);dy/dt=2xl 若找不到解析解，则返回其积分形式。若找不到解析解，则返回其积分形式。l 微分方程中用微分方程中用 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在在MatlabMatlab中的命令窗口中输入下面的命令中的命令窗口中输入下面的命令 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函数求解一些微分方程的通解