微分方程数值解课程设计教学文稿

微分方程数值解课程

设计

课程设计说明书

题目：常微分方程数值解法课程设计

学院（系）：理学院

年级专业：计算科学11-1

学生姓名：

指导教师：

教师职称：教授

燕山大学课程设计（论文）任务书

院（系）：理学院教学单位：仅供学习与交流，如有侵权请联系网站删除谢谢3

仅供学习与交流，如有侵权请联系网站删除谢谢4

年月日

燕山大学课程设计评审意见表

仅供学习与交流，如有侵权请联系网站删除谢谢5

仅供学习与交流，如有侵权请联系网站删除谢谢6

摘要

本文对常微分方程初值问题现有的数值解法进行了综述研究。主要讨论了几种常用的数值解法：即欧拉法，改进欧拉法，龙格库塔方法，阿达姆斯PECE,PMECME格式等。文章最后结合常见数值解法，对较为典型的微分方程模型进行数值求解，探讨了上述数值算法在实际建模问题中的应用。

论文阐述的是常微分方程数值解法的几个问题，通过对以下问题的求解

一.比较Adams四阶PECE模式和PMECME模式。

二.求解贝塞尔方程并与精确解比较。

三.小型火箭初始重量为1400kg，其中包括1080kg 燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N 的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。

来加强对用数值解法解常微分方程实际问题的能力。

关键词：常微分方程数值解 MATLAB

仅供学习与交流，如有侵权请联系网站删除谢谢i

目录

摘要 (i)

绪论 (1)

问题解答 (2)

总结 (18)

参考文献资料 (19)

仅供学习与交流，如有侵权请联系网站删除谢谢i

绪论

很多科学技术和工程问题常用微分方程的形式建立数学模型，因此微分方程的求解是很有意义的。建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些典型的方程，而对于绝大多数的微分方程问题，很难或者根本不可能得到它的解析解，实际问题终归结出来的微分方程主要靠数值解法。因此，研究微分方程求解的数值方法是非常有意义的。

仅供学习与交流，如有侵权请联系网站删除谢谢1

仅供学习与交流，如有侵权请联系网站删除 谢谢2

问题解答

问题一：比较Adams 四阶PECE 模式和PMECME 模式。

问题引出：将阿达姆斯方法显式与隐式方法作一对比，以说明预测——

把阿达姆斯外插公式与内插公式结合起来使用，先由前者提供初值，再由后者进行修正，所以Adams 预测-校正格式既利用了隐式方法较好的稳定性及精确性，有利用了显式方法的简易性，正是把两者结合起来，做到取长补短。 Adams 四阶预估-校正(PECE)公式：

(0)(1)(1)(1)(1)(1)

43321(0)(0)

444()(1)()(1)(1)(1)

43()()444[5559379(,u )C /24[195E (,u )

n n n n n n n n n n n n n n n n n n P u f f f f E f f t u h f f f f f f t ++++++++++++++++=+-+-??=??=++-+?

?=?１０４３２＋１１１ｕh/24］ ｕ９］ 而PMECME 模式公式为

(0)(1)(1)(1)(1)(1)

43321(0)(0)(1)(0)

4433(0)(0)

444()(1)()(1)(1)43'[55'59'37'9''[u ]'(,u')'/24['19'5'n n n n n n n n n n n n n n n n n n u f f f f M u u E f f t C u h f f f M E +++++++++++++++++=+-+-??=+-?

?=?=++-+?

???１０４３２P ｕh/24］ｕ251/270 ｕ９(1)(1)(1)(1)(0)4444()()444''[u ]'(,u')

n n n n n n n n f u u f f t +++++++=--=＋１１１］ｕ19/270 对于初值问题u ’=u-2t/u,u(0)=1,在单区间【0,1】上，用Adams 四阶预估-校正算法的PECE 模式及PMECME 模式求其数值解，取步长h=0.1利用计算结果估计数值解的局部误差主项。（真解u （t ）=sqrt(1+2t)）

编写matlab 程序：

function [ Un,e ] = pece()

%[ Un,e ] = pece()

f0=u-2*t/u;

v=[t,u];

U=zeros(1,11);

T=0:0.1:1;

f=zeros(1,11);

h=0.1;

U(1)=1;

f(1)=1;

for k=1:3

K1=subs(f0,v,[(k-1)*h,U(k)]);

K2=subs(f0,v,[(k-1)*h+1/2*h,U(k)+1/2*h*K1]);

K3=subs(f0,v,[(k-1)*h+1/2*h,U(k)+1/2*h*K2]);

K4=subs(f0,v,[(k-1)*h+h,U(k)+h*K3]);

U(k+1)=U(k)+1/6*h*(K1+2*K2+2*K3+K4);

f(k+1)=U(k+1)-2*T(k+1)/U(k+1);

end

for k=5:11

U(k)=U(k-1)+h/24*(55*f(k-1)-59*f(k-2)+37*f(k-3)-9*f(k-4)); f(k)=U(k)-2*T(k)/U(k);

U(k)=U(k-1)+h/24*(9*f(k)+19*f(k-1)-5*f(k-2)+f(k-3));

f(k)=U(k)-2*T(k)/U(k);

end

Un=U;

for k=1:11

zz(k)=sqrt(1+2*T(k));

end

e=U-zz;

end

function [ Un,e ] = pmecme()

%[ Un,e ] = pmecme()

syms t u

f0=u-2*t/u;

v=[t,u];

U=zeros(1,11);

T=0:0.1:1;

f=zeros(1,11);

h=0.1;

U(1)=1;

仅供学习与交流，如有侵权请联系网站删除谢谢3

f(1)=1;

for k=1:3

K1=subs(f0,v,[(k-1)*h,U(k)]);

K2=subs(f0,v,[(k-1)*h+1/2*h,U(k)+1/2*h*K1]);

K3=subs(f0,v,[(k-1)*h+1/2*h,U(k)+1/2*h*K2]);

K4=subs(f0,v,[(k-1)*h+h,U(k)+h*K3]);

U(k+1)=U(k)+1/6*h*(K1+2*K2+2*K3+K4);

f(k+1)=U(k+1)-2*T(k+1)/U(k+1);

end

for k=5:11

U(k)=U(k-1)+h/24*(55*f(k-1)-59*f(k-2)+37*f(k-3)-9*f(k-4)); U(k)=U(k)+251/270*(U(k)-U(k-1));

f(k)=U(k)-2*T(k)/U(k);

U(k)=U(k-1)+h/24*(9*f(k)+19*f(k-1)-5*f(k-2)+f(k-3));

U(k)=U(k)-19/270*(U(k)-U(k-1));

f(k)=U(k)-2*T(k)/U(k);

end

Un=U;

for k=1:11

zz(k)=sqrt(1+2*T(k));

end

e=U-zz;

end

运行结果及图形显示：

运行得：

>> [Un,e]=pece()

Un =

Columns 1 through 8

1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492

Columns 9 through 11

1.6125 1.6733 1.7321

e =

1.0e-05 *

Columns 1 through 8

0 0.0417 0.0789 0.1164 0.0571 0.0271 0.0127 0.0042

仅供学习与交流，如有侵权请联系网站删除谢谢4

Columns 9 through 11

-0.0013 -0.0054 -0.0088

pmecme

ans =

Columns 1 through 8

1.0000 1.0954 1.1832 1.2649 1.3398 1.4104 1.4773 1.5409

Columns 9 through 11

1.6016 1.6595 1.7147

编写a.m用于画图显示计算结果，

t0=0:0.001:1;

t1=0:0.1:1;

U0=sqrt(1+2*t0);

U2=pece();

U1=pmecme();

hold on

plot(t0,U0,'k')

plot(t1,U1,'--r','Marker','x')

plot(t1,U2,'--g','Marker','x')

legend('U=sqrt(1+2*t)','U1',’U2’'Location','NorthWest')

grid on

xlable('t')

ylable('Un')

hold off

运行结果如下：

仅供学习与交流，如有侵权请联系网站删除谢谢5

仅供学习与交流，如有侵权请联系网站删除

谢谢6

放大如下图

可知PMECME 模式比PECE 模式更为精确

问题二：求解贝塞尔方程并与精确解比较。

求解 x 2y ’’+xy ’+(x 2-n 2)y=0