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

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

《二阶常微分方程的数值求解.ppt》由会员分享，可在线阅读，更多相关《二阶常微分方程的数值求解.ppt（19页珍藏版）》请在冰点文库上搜索。

二阶常微分方程的数值求解,一.教学要求,掌握利用降阶把二阶常微分方程转化为一阶微分方程组，再利用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;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.m,RK_sys1.m,L2=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求解如下初值问题,解：

当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）;%K1andL1,rightf_sys1.m,RK_sys2.m,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）;%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：

求微分方程的通解，并验证。

y=dsolve（Dy+2*x*y=x*exp（-x2）,x）,symsx;diff（y）+2*x*y-x*exp（-x2）,利用dsolve函数求微分方程解析解,几点说明,如果省略初值条件，则表示求通解；,如果省略自变量，则默认自变量为t,dsolve（Dy=2*x,x）;dy/dx=2xdsolve（Dy=2*x）;dy/dt=2x,若找不到解析解，则返回其积分形式。

微分方程中用D表示对自变量的导数，如：

Dyy；D2yy；D3yy,例5：

求微分方程在初值条件下的特解，并画出解函数的图形。

y=dsolve（x*Dy+y-exp（x）=0,y

（1）=2*exp

（1）,x）ezplot（y）;,例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函数求解一些微分方程的通解,