常微分方程的数值解.doc

实验4 常微分方程的数值解 几道解方程的题——重要 1. x2y’’+xy’+(x2-n2)y=0 y=2,y’=- (贝赛尔方程，令n=0.5)，精确解y=sin x 解：首先将高阶方程装降阶化简为一阶常微分方程组 令y1=y’，令y2=y1’=y 将n=0.5代入，则原方程转化为： y1’=y2 y1=2, y2=- function dy=fangcheng(x,y) dy=[y(2);-y(2)/x-(1-1/(4*x^2))*y(1)]; xs=pi/2:0.1:100; y0=[2,-2/pi]; [x,y]=ode45(@fangcheng,xs,y0); [x,y] plot(x,y(:,1)),grid 2001 2． 设 用数值解法算出 y(1)= （精确到4位小数）, 你用的方法是 ，调用的 Matlab命令是 ，算法精度为 。 function dy=fangcheng2(x,y) dy=[y(2);y(2)*sin(x)-y(1)*exp(x)]; xs=0:0.01:1; y0=[1,0]; [x,y]=ode45(@fangcheng2,xs,y0); [x,y] 2000 1. 设用数值解法算出 y(1)= ,你用的方法是 ，调用的 Matlab命令是 ，算法精度为 。 function dy=fangcheng1(x,y) dy=[y(2);y(1)*sin(x)]; xs=[0,1]; y0=[1,0]; [x,y]=ode45(@fangcheng1,xs,y0); [x,y] 题目1 1. x2y’’+xy’+(x2-n2)y=0 y=2,y’=- (贝赛尔方程，令n=0.5)，精确解y=sin x 解：首先将高阶方程装降阶化简为一阶常微分方程组 令y1=y’，令y2=y1’=y 将n=0.5代入，则原方程转化为： y1’=y2 y1=2, y2=- （1）用向前欧拉公式： y1(n+1)=y1(n)+h*y2(n) y2(n+1)=y2(n)+h*[] y1(0) =2，y2(0)= -，x(n+1)= +n*h，h为步长 编程如下： clear all x=[pi/2:0.1:pi/2+100-0.05]; y1=2; y2=-2/pi; for i=1:999 y1(i+1)=y1(i)+0.1*y2(i); y2(i+1)=y2(i)-0.1*(y2(i)/x(i)+(1-0.25/x(i)^2)*y1(i)); end plot(x,y1),grid 但如果步长选得不一样(横坐标都从开始，到100左右结束，但中间选的点也不太一样)，即使采用同样的迭代公式，绘出的曲线却有很大不同： 程序： clear all x=[pi/2:0.01*pi:30*pi]; y1=2; y2=-2/pi; for i=1:2950 y1(i+1)=y1(i)+0.01*pi*y2(i); y2(i+1)=y2(i)-0.01*pi*(y2(i)/x(i)+(1-0.25/x(i)^2)*y1(i)); end plot(x,y1),grid 我个人认为，可能是因为欧拉公式每算一次都会产生误差，如果取的点正好位置不太合适，会导致误差一步步增大，累加起来，就像蝴蝶效应一样，会产生和真实状况截然不同的结果。这也是用数值方法求解方程最怕出现的问题，也是应该解决的问题。这说明向前欧拉公式并不是一种很好的计算方法，误差较大（只有1阶精度），而且容易失真。这一点在下面还要说明。 （2）龙格—库塔方法 y1’=y2 y1=2, y2=- 编写如下M文件： function dx=suan(t,x) %建立名为suan的函数M文件 dx=[x(2);-x(2)/t-(1-0.25/t^2)*x(1)]; %表示方程组 用4级5阶龙格—库塔公式进行计算。编程如下： clear all ts=[pi/2:0.1:100]; %使用龙格—库塔方法 x0=[2,-2/pi]; [t,x]=ode45(@suan,ts,x0) plot(t,x(:,1)),grid,title(龙格-库塔方法),pause %绘出精确解y=sin x图像 y=sin(t).*sqrt(2*pi./t) plot(t,y),grid,title(精确解) 绘出曲线如下： 比较分析：从曲线上可以看出，在