数值分析课程设计常微分方程组初值问题数值解的实现和算法分析毕业论文Word文档下载推荐.docx
《数值分析课程设计常微分方程组初值问题数值解的实现和算法分析毕业论文Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《数值分析课程设计常微分方程组初值问题数值解的实现和算法分析毕业论文Word文档下载推荐.docx（15页珍藏版）》请在冰豆网上搜索。
（2）
3.2.1用Runge-Kutta方法计算解决一阶微分方程组初值问题的基本思路
（2）式形式上与常数微分方程初值问题是一样的，只要注意向量函数运算与其表示，就可以用初值问题的求解格式得到常微分方程组初值问题
（2）的求解格式，由初值问题的经典Runge-kutta公式可得一阶常微分方程组初值问题
（2）的Runge-kutta公式：
注意上式是向量形式，其对应的分量形式为：
微分方程理论告诉我们，高阶微分方程可转化为一阶微分方程组来研究，因此可以用一阶微分方程组初值问题揭发来解高阶微分方程初值问题。
高阶微分方程初值问题的形式为：
（3）
令则
（2）化为了一阶微分方程组初值问题：
Runge-kutta方法巧妙利用函数
在一些点上的函数值的线性组合，获得了高阶的数值解法，它避开了要获得高阶方法须对
求高阶导数的不便，是离散化方法中Tayl情况，其中在准确性的工作量的综合效果看，经典的Runge-kutta方法是首选or展开法的一个应用。
Runge-kutta方法主要用于定步长的。
Runge-kutta方法也常用于对多步法提供初值。
3.2.2用改进Euler方法计算解决一阶微分方程组初值问题的基本思路
改进Euler方法需要用Euler方法求出一个预测值
然后再用梯形公式校正一次得到
，即所求结果的迭代格式。
（4）
为了方便编程可将（4）式改变为如下格式
4用matlab语言编程解决相关问题
4.1四阶Runge-kutta方法的Matlab编程实现
function[T]=Runge_Kutta（f,x0,y0,h,n）
%
[T]=Runge_Kutta（f,x0,y0,h,n）
f
待解方程（组）
x0
初试自变量值
y0
初试函数值
h
步长
n
步数
ifnargin<
5
n=100;
end
r=size（y0）;
r=r
（1）;
s=size（x0）;
s=s
（1）;
r=r+s;
T=zeros（r,n+1）;
T（:
1）=[y0;
x0];
fort=2:
n+1
k1=feval（f,[T（1:
r-s,t-1）;
x0]）;
k2=feval（f,[k1*（h/2）+T（1:
x0+h/2]）;
k3=feval（f,[k2*（h/2）+T（1:
k4=feval（f,[k3*h+T（1:
x0+h]）;
x0=x0+h;
t）=[T（1:
r-s,t-1）+（k1+k2*2+k3*2+k4）*（h/6）;
End[2]
4.2Euler改进方法Matlab编程实现
用Euler改进方法编写Matlab程序为：
function[x,y]=eulerpro（fun,x0,xfinal,y0,n）;
h=（xfinal-x0）/n;
x
（1）=x0;
y
（1）=y0;
fori=1:
n
x（i+1）=x（i）+h;
y1=y（i）+h*feval（fun,x（i）,y（i））;
y2=y（i）+h*feval（fun,x（i+1）,y1）;
y（i+1）=（y1+y2）/2;
5编程解决
5.1输入计算题目
functiondY=dydx（X,Y）
dY
（1）=-0.013*Y
（1）-1000*Y
（1）*Y
（2）;
dY
（2）=-2500*Y
（2）*Y（3）;
dY（3）=0.013*Y
（1）-1000*Y
（1）*Y
（2）-2500*Y
（2）*Y（3）
dY=[dY
（1）;
dY
（2）;
dY（3）];
5.2用Runge-Kutta方法的Matlab编程解法
function[k,X,Y,wucha,P]=RK4z（dydx,a,b,CT,h）
n=fix（（b-a）/h）;
X=zeros（n+1,1）;
Y=zeros（n+1,length（CT））;
X=a:
h:
b;
Y（1,:
）=CT'
;
fork=1:
k1=feval（dydx,X（k）,Y（k,:
））
x2=X（k）+h/2;
y2=Y（k,:
）'
+k1*h/2;
k2=feval（dydx,x2,y2）;
k3=feval（dydx,x2,Y（k,:
+k2*h/2）;
k4=feval（dydx,X（k）+h,Y（k,:
+k3*h）;
Y（k+1,:
）=Y（k,:
）+h*（k1'
+2*k2'
+2*k3'
+k4'
）/6;
k=k+1;
fork=2:
wucha（k）=norm（Y（k）-Y（k-1））;
k=k+1;
X=X（1:
n+1）;
Y=Y（1:
n+1,:
）;
k=1:
n+1;
wucha=wucha（1:
k,:
P=[k'
X'
Y,wucha'
];
调用dydx.m求解，在MATLAB工作窗口输入程序
CT=[1;
1;
0];
h=0.0001;
[k,X,Y,wucha,P]=RK4z（dydx,0,0.01,CT,h）,
H=[0.0001,0.01];
[x,y]=ode15s（'
dydx'
H,CT）;
plot（X,Y（:
1）,'
g-'
x,y（:
bo'
X,Y（:
2）,'
m:
'
cp'
3）,'
r.'
kd'
）
xlabel（'
轴\itx'
ylabel（'
轴\ity'
title（'
分别用自定义函数和ode15s函数求解刚性方程方程组的图形'
legend（'
用RK4z函数解刚性方程的y1的曲线'
'
用ode15s函数解刚性方程的y1的曲线'
用RK4z函数解刚性方程的y2的曲线'
用ode15s函数解刚性方程的y2的曲线'
'
用RK4z函数解刚性方程的y3的曲线'
用ode15s函数解刚性方程的y3的曲线'
5.3用改进Euler方法的Matlab编程解法
function[x,y]=eulerpro（dydx,a,b,CT,h）;
x=zeros（n+1,1）;
y=zeros（n+1,length（CT））;
y（1,:
f1=feval（dydx,x（i）,y（i,:
））;
y1（i+1,:
）=y（i,:
）+h*f1'
f2=feval（dydx,x（i+1）,y1（i+1,:
y2（i+1,:
）+h*f2'
y（i+1,:
）=（y1（i+1,:
）+y2（i+1,:
））/2;
a=0;
b=0.01;
[X,Y]=eulerpro（dydx,a,b,CT,h）,
dy123'
用eulerpro函数解刚性方程的y1的曲线'
用eulerpro函数解刚性方程的y2的曲线'
用eulerpro函数解刚性方程的y3的曲线'
6计算结果
运行结果详见附录
6.1用四阶Runge-Kutta方法的Matlab编程解法的结果以与与精确解的比较
6.2用改进Euler方法的Matlab编程解法的结果以与与精确解的比较
（注：
计算结果详见附录1，附录2）
7.结果分析
由图可知次方法与精确解的契合度非常好，基本上与精确解保持一致，由此可见四阶Runge-Kutta方法是一种高精度的单步方法。
此法对比改进Euler方法精确度更高。
但是，相对的计算步骤比Euler改进方法要繁琐。
综上，当计算低精度问题时可以使用改进Euler方法来处理问题，而如果精度要求较高，就要使用四阶Runge-Kutta方法。
致
本次课程设计主要针对一阶常微分方程组的初值问题，利用改进Euler方法以与Runge-Kutta方法解决一阶产微分方程的算法进行推广到一阶常微分方程组的算法。
其中一部分编程思路参考了相关的参考书，经过反复验证然后得到的结果。
其中在编写改进Euler方法的程序时，由于参考程序是非方程组的程序，在改写的时候忽略了向量的方向导致不匹配，通过仔细检查后成功运行程序。
回顾起此次课程设计，至今我仍感慨颇多，的确，从找参考，设计到定稿，从理论到实践，在整整一星期的日子里，可以说得是苦多于甜，但是可以学到很多很多的的东西，同时不仅可以巩固以前所学过的知识，而且学到了很多在书本上所没有学到过的知识。
通过这次课程设计使我懂得了理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，从理论中得出结论，才能提高自己的实际动手能力和独立思考的能力。
同时在设计的过程中发现了自己的不足之处，对以前所学过的知识理解得不够深刻，掌握得不够牢固。
通过这次课程设计之后，一定把以前所学过的知识仔细复习，特别是要注意细节部分。
参考文献
[1]庆扬,王能超，易大义.数值分析[M].清华大学,2008.
[2]宋叶志，贾东永.Matlab数值分析与应用[M].机械工业,2009.
[3]黄明游,果忱.数值分析[M].高等教育,2008.
附录
改进Euler方法计算结果
Y=
1.00001.00000
0.90501.0125-0.0825
0.81721.0427-0.1401
0.73511.0865-0.1785
0.65791.1405-0.2016
0.58551.2020-0.2125
0.51771.2686-0.2137
0.45471.3378-0.2075
0.39651.4075-0.1959
0.34351.4757-0.1808
0.29561.5408-0.1636
0.25271.6016-0.1456
0.21491.6573-0.1278
0.18171.7074-0.1108
0.15301.7518-0.0952
0.12831.7906-0.0811
0.10721.8243-0.0686
0.08931.8531-0.0576
0.07421.8776-0.0482
0.06151.8983-0.0402
0.05091.9157-0.0334
0.04201.9302-0.0277
0.03471.9424-0.0229
0.02861.9525-0.0189
0.02351.9608-0.0156
0.01941.9678-0.0129
0.01591.9735-0.0106
0.01311.9782-0.0087
0.01081.9821-0.0072
0.00881.9853-0.0059
0.00731.9879-0.0048
0.00601.9901-0.0040
0.00491.9918-0.0033
0.00401.9933-0.0027
0.00331.9945-0.0022
0.00271.9955-0.0018
0.00221.9963-0.0015
0.00181.9970-0.0012
0.00151.9975-0.0010
0.00121.9979-0.0008
0.00101.9983-0.0007
0.00081.9986-0.0005
0.00071.9989-0.0004
0.00061.9991-0.0004
0.00051.9992-0.0003
0.00041.9994-0.0002
0.00031.9995-0.0002
0.00021.9996-0.0002
0.00021.9996-0.0001
0.00021.9997-0.0001
0.00011.9998-0.0001
0.00011.9999-0.0001
0.00011.9999-0.0000
0.00001.9999-0.0000
0.00002.0000-0.0000
………
0.00002.0000-0.0000
四阶Runge-Kutta方法结算结果
P=
1.000001.00001.000000
2.00000.00010.90451.0112-0.08430.0955
3.00000.00020.81641.0401-0.14350.0881
4.00000.00030.73431.0814-0.18430.0821
5.00000.00040.65751.1311-0.21140.0769
6.00000.00050.58551.1863-0.22820.0719
7.00000.00060.51851.2445-0.23690.0670
8.00000.00070.45651.3042-0.23930.0620
9.00000.00080.39951.3638-0.23670.0570
10.00000.00090.34751.4222-0.23020.0520
11.00000.00100.30061.4786-0.22070.0469
12.00000.00110.25861.5324-0.20900.0420
13.00000.00120.22131.5830-0.19570.0373
14.00000.00130.18841.6302-0.18140.0328
15.00000.00140.15971.6737-0.16660.0287
16.00000.00150.13481.7134-0.15170.0249
17.00000.00160.11341.7495-0.13700.0214
18.00000.00170.09501.7820-0.12290.0184
19.00000.00180.07941.8111-0.10950.0156
20.00000.00190.06621.8368-0.09700.0132
21.00000.00200.05501.8596-0.08540.0112
22.00000.00210.04561.8796-0.07470.0094
23.00000.00220.03781.8971-0.06510.0079
24.00000.00230.03121.9123-0.05650.0066
25.00000.00240.02581.9254-0.04880.0055
26.00000.00250.02121.9367-0.04200.0045
27.00000.00260.01751.9465-0.03600.0037
28.00000.00270.01441.9548-0.03080.0031
29.00000.00280.01181.9619-0.02620.0026
30.00000.00290.00971.9680-0.02230.0021
31.00000.00300.00801.9731-0.01890.0017
32.00000.00310.00661.9775-0.01600.0014
33.00000.00320.00541.9811-0.01350.0012
34.00000.00330.00441.9842-0.01130.0010
35.00000.00340.00361.9868-0.00950.0008
36.00000.00350.00301.9890-0.00800.0007
37.00000.00360.00241.9908-0.00670.0005
38.00000.00370.00201.9924-0.00560.0004
39.00000.00380.00161.9937-0.00470.0004
40.00000.00390.00131.9947-0.00390.0003
41.00000.00400.00111.9956-0.00330.0002
42.00000.00410.00091.9964-0.00270.0002
43.00000.00420.00071.9970-0.00230.0002
44.00000.00430.00061.9975-0.00190.0001
45.00000.00440.00051.9979-0.00160.0001
46.00000.00450.00041.9983-0.00130.0001
47.00000.00460.00031.9986-0.00110.0001
48.00000.00470.00031.9988-0.00090.0001
49.00000.00480.00021.9990-0.00070.0000
50.00000.00490.00021.9992-0.00060.0000
51.00000.00500.00011.9993-0.00050.0000
52.00000.00510.00011.9994-0.00040.0000
53.00000.00520.00011.9995-0.00030.0000
54.00000.00530.00011.9996-0.00030.0000
55.00000.00540.00011.9997-0.00020.0000
56.00000.00550.00011.9997-0.00020.0000
57.00000.00560.00001.9998-0.00020.0000
58.00000.00570.00001.9998-0.00010.0000
59.00000.00580.00001.9998-0.00010.0000
60.00000.00590.00001.9999-0.00010.0000
61.00000.00600.00001.9999-0.00010.0000
62.00000.00610.00001.9999-0.00010.0000
63.00000.00620.00001.9999-0.00000.0000
64.00000.00630.00001.9999-0.00000.0000
65.00000.00640.00001.9999-0.00000.0000
66.00000.00650.00001.9999-0.00000.0000
67.00000.00660.00002.0000-0.00000.0000
………………
101.00000.01000.00002.0000-0.00000.0000
翻译
关于feval的Matlabhelp原文如下：
FEVALExecutethespecifiedfunction.
FEVAL（F,x1,...,xn）evaluatesthefunctionspecifiedbyafunction
handleorfunctionname,
升级会员 