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

二阶常微分方程的数值求解一.教学要求掌握利用降阶把二阶常微分方程转化为一阶微分方程组，再利用Euler方法数值求解,并能利用MATLAB软件进行数值计算和符号运算。二.教学过程考虑如下的二阶微分方程初值问题若令'zy,则上述初值问题可以转化为如下一阶微分方程组初值问题2012(,,'),(),'(),[,]dyfxyyyayyayxabdx100'()(),[,]'()(,(),()),[,](),()yxzxxabzxfxyxzxxabzayzyay利用Euler方法求解上述方程组可得如下数值格式0011112(),'(),(,,),,.kkkkkkkkkkyayyazyyhzzzhfxyzkxxh()'()kkkkyyxzyx其中是的近似，是的近似利用四阶R-K方法求解上述方程组可得如下数值格式11234112341121211323224343322622622222222(),(),,(,,),,(,,),,(,,),,(,,).kkkkkkkkkkkkkkkkkkkkhyyKKKKhzzLLLLKzLfxyzhhhhKzLLfxyKzLhhhhKzLLfxyKzLkKzhLLfxhyhKzhL12,()'()kkkkyyxzyx其中是的近似，是的近似例1：用Euler法求解如下初值问题22020101[,](),'()dyyxdxyy当h=0.1，即n=20时，Matlab源程序见Euler_sys1.m解：该问题的真解为xye002020101','()(),[,]'()(),[,](),().令则该初值问题可以转化为zyyxzxxzxyxxzzyclc;clear;h=0.1;a=0;b=2;x=a:h:b;y(1)=1;z(1)=-1;fori=1:length(x)-1y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*y(i);endplot(x,y,'r+',x,exp(-x),'k-');xlabel('Variablex');ylabel('Variabley');Euler_sys1.m数值解与真解如下图例2：利用4阶R-K方法求解例1，并与Euler方法进行比较。解当h=0.1，即n=20时，R-K方法的Matlab源程序见RK_sys1.m，数值结果见下图functionw=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;%初值fori=1:length(x)-1%****EulerMethod****%Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);Euler_z(i+1)=Euler_z(i)+h*Euler_y(i);%*****R-K4Method*****%K1=RK_z(i);L1=rightf_sys1(x(i),RK_y(i),RK_z(i));%K1andL1K2=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);%K2andL2K3=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);%K3andL3K4=RK_z(i)+h*L3;L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);%K4andL4RK_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('Variablex');ylabel('Variabley');例3：分别用Euler法和R-K4求解如下初值问题2221020101()[,](),'()xdyyexxdxyy解：该问题的真解为cosxyxxe00221020101','()(),[,]'()()(),[,](),().令则该初值问题可以转化为xzyyxzxxzxyxexxzzy当h=0.1，即n=20时，Matlab源程序见RK_sys2.m,数值结果如下图functionw=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;fori=1:length(x)-1%****EulerMethod****%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-K4Method*****%K1=RK_z(i);L1=rightf_sys2(x(i),RK_y(i),RK_z(i));%K1andL1rightf_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);%K2andL2K3=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);%K3andL3K4=RK_z(i)+h*L3;L4=rightf_sys2(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);%K4andL4RK_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('Variablex');ylabel('Variabley');dsolve的调用格式y=dsolve('eq1','eq2',...,'cond1','cond2',...,'v')其中y为输出，eq1、eq2、...为微分方程，cond1、cond2、...为初值条件，v为自变量，如果不指定v作为自变量，则默认t为自变量。例4：求微分方程的通解，并验证。22xdyxyxedxy=dsolve('Dy+2*x*y=x*exp(-x^2)','x')symsx;diff(y)+2*x*y-x*exp(-x^2)利用dsolve函数求微分方程解析解几点说明如果省略初值条件，则表示求通解；如果省略自变量，则默认自变量为tdsolve('Dy=2*x','x');％dy/dx=2xdsolve('Dy=2*x');％dy/dt=2x若找不到解析解，则返回其积分形式。微分方程中用D表示对自变量的导数，如：Dyy'；D2yy''；D3yy'''例5：求微分方程在初值条件下的特解，并画出解函数的图形。0'xxyyey=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')ezplot(y);12()ye2220100cos()(),'()求二阶常微分方程的通解dyxydxyy例6在Matlab中的命令窗口中输入下面的命令symsxyS=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函数求解一些微分方程的通解