数学实验作业四

1、数学实验作业四题目：P120. 1.d); 4); 5)日期:2003-3-23【实验内容】：用欧拉方法和龙格-库塔方法求微分方程数值解,画曲解的图形，对结果进行比较分析。以下方程供选择:a)y' = y +2x, y(0) =1,0 Ex E1,精确解 y =3ex2x 2.b)22y =x -y ,y(0) =0,或y(0) =1.c)x2y''+xy' + (x2n2)y=°,yC)=2,y'(2) = q(Bessel 方程，令n=1/2),精确解2 二 y=sinx x3d) y'' y cosx =0,y(0) =1

2、,y'(0) =0,哥级数解,1 22 49 6 55 8y = 1 - x '' x - x "x2!4!6!8!【模型建立】：原方程化为 y'' = -ycosx,y(0) fy'B)=。x(l) = y记凶2)=y',并用t代替x,则原方程化为:x(1)' =x(2)x(1)建=1x(2)' = x(1)cost .日x(2)|t包=0 .于是可以用龙格-库塔法求解。【模型求解】：用Matlab作龙格库塔法求解：% chapter 4 -1.d%此函数是微分方程组function Xdot=ch41dfun

3、(t,x)Xdot=x(2),-cos*x'%ch41d.mfunction I=ch41d(a)x0=1,0't,x=ode45('ch41dfun',0,a,x0);y=x(:,1);plot(t,y,'r');hold on;gtext('y');y1=1-1/jiecheng(2)*t.A2+2/jiecheng(4)*t.A4-9/jiecheng(6)*t.A6+55/jiecheng(8)*t.A8;plot(t,y1,'b');gtext('y1');t,y,yihold off;运

4、行程序可以得到：(只保留图像，数据省略)图1>> ch41d(10)(数据省略，yl图象也省略)图25>> ch41d(15)>> ch41d(20)图47【结果分析】由图1得：y（龙格库塔方法）和y1（级数近似解）在。到 大约1.5的区间内是完全吻合的，从 x=1.5之后，两条曲线 开始分离。之后y的变化趋势可见图2-图4,呈幅度越来 越大的上下震荡。对y取绝对值后再取对数，得到其图象为：由此图可见：|y|基本是随x指数增加的，并上下振动，可 以猜测|y|接近于ex+osx的形式。四：求解“火箭发射”问题。【模型分析】从上图可以看生，开始的几个点基本呈线性

5、； 而以后的几模型假设:1) xwx0时，F 与 x 呈正比，F(x) = kx;模型建立:采用最小二乘法，目标函数为模型求解:为求生J的最小值，应保证 编程语句如下：function k,a,b,c=number8(m) x = 1,2, 4, 7, 9,12,13,15,17;运行程序可以得到:x0 = 7.9342同时得到拟合曲线如下:结论:1.7000x x < 7,9342F =2-0,097仅2 3.2404x - 6.0417 x _ 7,9342五：一只小船渡过宽为 d的河流（见下图），目标是起点 A正对着的另一岸 B点。已知河水流速 v1与船在静水中的速度v2之比为ko

6、a)建立小船航线的方程，求其解析解。b)设d=100m,v1=1m/s,v2=2m/s,用数值解法求渡河所需时问、任意时刻小船的位置及航行曲线，作图，并与解析解比较。c)若流速v1为0,0.5,1.5,2(m/s),结果将如何？【模型建立】：如图，以B为原点，沿河岸向右为 x轴正向，垂直河岸 向下为y轴正向，建立坐标系。设在 t时刻，船在x方向上 的位移是x(t),在y方向上的位移是y(t),则：在t时刻，船在 X方向上的速度是X，在y方向上的速度是y'，将船的速 度V和静水速度V1在x,y轴方向上分解，可得：Vx =V1 -V2sin FVy - -V 2cos 1又因为船头始终指向

7、Bxtan =一点，所以 y,所以:Vx =V1 -V2Vy = -V2xx2 y2yx x'(t) =V1 -V2 x2 y2 y'(t) =N2= xy2这就是本题的微分方程。初始条件为:x(0)=0,y(0)=-d【模型求解】：1.解析解：x ;：；cosu令1y = Psine ,将直角坐标化成极坐标,1 P'cos - Psin 二二'=V1-V2cos1I i 'sin 二 、cosi 二'-V2sin1;则原微分方程化成化简后为:'=V1cos -V2H'； '=V1sinu消去t,得:'V2二(-c

8、tg V1于是解得:一2sin日，其中V2R =V12.数值解法：根据上述微分方程，用龙格-库塔方法求解:程序如下：%ch45fun.m%第四章第5题%此函数是微分方程组function Xdot=ch45fun(t,x,V1,V2)%x(1)代表x;x(2)代表y%加入限制条件，防止无限循环if (norm(x)>1e-5)Xdot=V1-V2*x(1)/sqrt(x(1)A2+x(2)A2),-V2*x(2)/sqrt(x(1)A2+x(2F2)'elseXdot=0,0'end%ch45.m%第四章第 5 题function I=ch45(V1,V2,d,a)x0=

9、0,-d'% 初始条件t,x=ode45(ch45fun,0,a,x0,V1,V2);plot(x(:,1),x(:,2),'r');hold on;%作图t,x(:,1),x(:,2)%打印数据%下面是作出精确解的图象seta=linspace(-pi/2,0,100);rou=d*(abs(tan(seta/2)(V2/V1)./sin(seta);xp=-rou.*cos(seta);yp=-rou.*sin(seta);plot(xp,yp,'b');hold off;需要说明：在ch45fun.m 中，加入了 (norm(x)>1e-5)

10、 的限制条件，以保证在船离B 点足够近时中止运算。否则无法得出正确结果。运行结果：(中间数据省略)可知渡河时航线如图，所用时间为： 66.7秒当v1=0,0.5,1.5,2时，结果分别为（图像）:v1=0 10 20 30 40-50 BO 70-80 901001-0.8-0 6 -a J-口N00 20.4 O S 0.81说明在静水中，船沿直线到达B点。这与直观经验相符合渡河时间为50秒；v1=0.5渡河时间为：53.3秒v1=1.517渡河时间为：114.3秒v1=2渡河时间：从图上看由，到t=1000秒时，船已到达对岸,但是并没有到达 B点，而是在B点下游50米处。由于船头 指向B点

11、，即船头指向逆流方向， 且船速（静水）等于水速, 可知船将保持原地不动。也就是说，船永远到达不了B点。所以渡河时间为无穷大。v1=2.5渡河时间：与v1=2时情况类似，船能到达对岸，但是是在 B点下游。由于船速（静水）小于水速，船将被水冲得顺流 而下，同样永远到达不了 B点。【结果分析】1 . 以上图中红线是龙格库塔方法的结果，蓝线是解析解的结果。从图中可以看出，两种方法的结果基本上是相符的。但在接近 B 点时，解析解将无法得到正确解（因为会出现分母为零的情况，而Matlab 不会计算极限） 。当 V1=0 时，解析解也得不到正确解（同样是因为出现了分母为零） 。2 .随着V1的变化，船的航线也在变化。当