下载文档原格式
课题10. 常微分方程初值问题的数值方法一．问题提出（1）利用欧拉方法和改进的欧拉方法求解初值问题：dy/dx=4*x/y-x*y; y(0)=3; 其中0<x<=2二.问题算法、c语言编程和上机运算结果1.欧拉法算法：输入微分方程f（x，y）；输入积分部数n；输入初值x0，y0；输入步长h；利用k1=f(xn,yn)y（n+1）=yn+k1*h；n=0,1,2……采取不断循环计算；输出x1，x2，……xn.程序：/*微分方程*/float f(x,y)float x,y;{float z;z=4*x/y-x*y;return(z);}/* 欧拉法*/float EULAR(f)float (*f)();{int i,n;float x0,y0,x,y,k1,h;printf("\n请输入积分步数n:");scanf("%d",&n);printf("\n请输入初值x(0) y(0):");scanf("%f%f",&x0,&y0);printf("\n请输入步长h:");scanf("%f",&h);printf("\n x y"); printf("\n %f %f",x0,y0); for(i=1;i<=n;i++){x=x0+h;k1=(*f)(x0,y0);y=y0+h*k1;printf("\n %f %f",x,y); x0=x;y0=y;}}main(){EULAR(f);}结果：(1)h=0.1时（2）h=0.2时（3）h=0.4时2.改进的欧拉法算法：输入微分方程f（x，y）；输入积分部数n；输入初值x0，y0；输入步长h；利用k1=f(xn,yn)k2=f(xn+h,yn+h*k1)y（n+1）=yn+（k1+k2）/2；n=0,1,2……采取不断循环计算；输出x1，x2，……xn.程序：/*微分方程*/float f(x,y)float x,y;{float z;z=y-2.0*x/y;return(z);}/* 改进欧拉法*/float EULAR(f)float (*f)();{int i,n;float x0,y0,x,y,k1,k2,h;printf("\n请输入积分步数n:");scanf("%d",&n);printf("\n请输入初值x(0) y(0):");scanf("%f%f",&x0,&y0);printf("\n请输入步长h:");scanf("%f",&h);printf("\n x y"); printf("\n %f %f",x0,y0); for(i=1;i<=n;i++){x=x0+h;k1=(*f)(x0,y0);k2=(*f)(x,y0+h*k1);y=y0+h*(k1+k2)/2.0;printf("\n %f %f",x,y); x0=x;y0=y;}}main(){EULAR(f);}结果：（1）当h=0.1时（2）当h=0.2时（3）当h=0.4时三.结果分析讨论1.对比欧拉法，改进的欧拉法和精确解的结果可知，改进的欧拉法所得到结果的精度比欧拉法的大，这是因为改进的欧拉法融入了属于隐式公式的梯形公式，它的计算数值解的精度要比欧拉公式好。
第九章常微分方程初值问题数值解法图9-1n 作为()n x y 的近似值，得 ()n n y x hf ,)y x ,两边从n x 到1+n x 积分，得()dx x y x f x y x n nx x n n ⎰+=-+1))(,()1 矩形公式计算上式右侧积分，即()()x x x x x d x y x f dx x y x f n nn n⎰⎰++≈11,))(,()n ，得()n n n n y x hf y y ,1+=+，故欧拉法也称为矩形法。
为了达到较高精度的计算公式，对欧拉法进行改进，用梯形公式计算()()([1,2))(,(1++≈+n n n x f x y x f hdx x y x f n 的近似值，得9.2 龙格—库塔法前面讨论的欧拉法与改进的欧拉法都是一步法，即计算y 1+n 时，只用到前一步值。
龙格—库塔(Runge-Kutta)法(简称为R-K 方法)不是通过求导数的方法构造近似公式，而是通过计算不同点上的函数值，并对这些函数值作线性组合，构造近似公式，再把近似公式与解的泰勒展开式进行比较，使前面的若干项相同，从而使近似公式达到一定的阶数。
我们先分析欧拉法与预估—校正法。
对于欧拉法⎩⎨⎧=+=+),(111n n n n y x hf k k y y 每步计算f 的值一次，其截断误差为O (2h )。
对于预估—校正法()()⎪⎪⎩⎪⎪⎨⎧++==++=+121211,,2121k y h x hf k y x hf k k k y y n n n n n n 每步计算f 的值两次，其截断误差为O (3h ).下面对预估—校正法进行改进，将该公式写成更一般的形式()()bh y ah x hf k y x hf k k R k R y y n n n n n n ++==++=+,,2122111 (2.1)其中b a R R ,,,21为待定常数。
选择这些常数的原则是在)(n n x y y =的前提下，使11)(++-n n y x y )的阶尽量高。
常微分方程的初值问题常微分方程是数学中的一种重要工具，它能够描述许多自然界和社会现象的变化规律。
而常微分方程的初值问题则是常微分方程研究中的常见问题之一，它需要确定未知函数及其导数在某个特定点的值。
本文将介绍常微分方程的初值问题的定义、求解方法以及实际应用。
一、初值问题的定义在常微分方程中，初值问题是指在已知微分方程的解的条件下，需要确定一个特定点上未知函数及其导数的值。
具体而言，考虑一个形如dy/dx=f(x,y)的一阶常微分方程，其中x是自变量，y是因变量，f是已知的函数。
若已知y(x0)=y0，则求解这个微分方程的过程即为解决初值问题。
二、求解方法对于常微分方程的初值问题，可以使用多种方法进行求解，下面将介绍两种常见的方法：欧拉方法和四阶龙格-库塔方法。
1. 欧拉方法欧拉方法是一种简单而直观的求解常微分方程的数值方法。
它的基本思想是将求解区间等分为多个小区间，然后通过逐步逼近的方式计算未知函数的近似值。
具体步骤如下：- 将求解区间[a, b]等分为n个小区间，步长h=(b-a)/n。
- 定义网格节点xi=a+i*h，i=0,1,2,...,n。
- 初始条件为y(x0)=y0，通过递推公式y(xi+1) = y(xi) + h*f(xi, y(xi))，计算出近似值y(xi+1)。
- 重复上述步骤，直到计算到需要的点。
欧拉方法的优点是简单易懂，但对于某些特定的微分方程，其数值解可能不够精确。
2. 四阶龙格-库塔方法四阶龙格-库塔方法是一种更为精确的求解常微分方程的数值方法，它通过计算多个逼近值的组合来提高计算精度。
具体步骤如下：- 将求解区间[a, b]等分为n个小区间，步长h=(b-a)/n。
- 定义网格节点xi=a+i*h，i=0,1,2,...,n。
- 初始条件为y(x0)=y0，通过递推公式计算逼近值k1、k2、k3和k4。
- k1 = h*f(xi, y(xi))- k2 = h*f(xi + h/2, y(xi) + k1/2)- k3 = h*f(xi + h/2, y(xi) + k2/2)- k4 = h*f(xi + h, y(xi) + k3)- 计算近似值y(xi+1) = y(xi) + (k1 + 2k2 + 2k3 + k4)/6。
贵州师范大学数学与计算机科学学院学生实验报告
课程名称： 数值分析 班级： 实验日期： 年 月 日
学 号： 姓名： 指导教师：
实验成绩：
一、实验名称
实验六: 常微分方程初值问题数值解法
二、实验目的及要求
1. 让学生掌握用Euler法, Runge-Kutta法求解常微分方程初值问题.
2. 培养Matlab编程与上机调试能力.
三、实验环境
每人一台计算机，要求安装Windows XP操作系统，Microsoft office2003、
MATLAB6.5(或7.0).
四、实验内容
1. 取步长h=0.1,0.05,0.01, ,用Euler法及经典4阶Runge-Kutta 法求解初值
问题
1)0()10(2222'y
tttyy
要求:
1) 画出准确解(准确解22teyt)的曲线,近似解折线;
2) 把节点0.1和0.5上的精确解与近似解比较,观察误差变化情况.
2. 用 Euler法,隐式Euler法和经典4阶R-K法取不同步长解初值问题
21)0(],1,0[,50'y
xyy
并画出曲线观察稳定性.
注：题1必须写实验报告
五、算法描述及实验步骤
Euler法：
输入 000),(,,,),,(yaxxhbayxf
输出 Euler解y
步1 ),,2,1(;mnhnaxhabmn
步2 对1,,2,1,0mn执行),(1nnnnyxfhyy
步3 输出
T
myyyy),,,(21
经典4阶R-K法：
输入 000),(,,,),,(yaxxhbayxf
输出 4阶R-K解y
步1 ),,2,1(;mnhnaxhabmn
步2 对1,,2,1,0mn执行),(1nnyxfK，)5.0,(15.02hKyxfKnn，
)5.0,(25.03hKyxfKnn，),(314hKyxfKnn
)22(643211KKKKhyynn
步3 输出
T
myyyy),,,(21
六、调试过程及实验结果
>> shiyan6
Y1 =
0.8000 0.6620 0.5776 0.5401 0.5441 0.5853 0.6602 0.7662 0.9009 1.0627
Y2 =
0.8287 0.7103 0.6388 0.6093 0.6179 0.6612 0.7366 0.8419 0.9753 1.1353
e1 = 0.0287 e2 = 4.2469e-006 e1 =
0.0738
e2 =
1.1609e-005
注：至于h=0.05、0.01的情况将程序中的h值作相应的改动即可得。
七、总结
在1阶的Euler法解初值问题时，若步长过大的话，误差将会较大，其解不可靠，只有
控制步长，在一定误差范围内才可用。
就精度阶而言经典4阶Runge-Kutta法最可取，但并不是阶高的方法就一定可用，还必
须考虑初值问题的性态、数值方法的计算量稳定性等因素。因此在实际计算中，应根据问题
的具体情况来选择合适的方法。
八、附录（源程序清单）
图像、计算、误差比较程序：
x=0:0.001:1; y=exp(-2*x)+x.^2; plot(x,y,'r') axis([0,1,0.5,1.2]) hold on a=0;b=1;y0=1;h=0.1;d=y0; Y1=Euler('fun',a,b,y0,h) u1=0:0.001:h; v1=Y1(1)+0*u1; plot(u1,v1,'g--') hold on u2=0:0.001:5*h; v2=Y1(5)+0*u2; plot(u2,v2,'g--') hold on Y1=[d,Y1]; t=0:h:1; scatter(t,Y1,'r') hold on
plot(t,Y1)
hold on
Y2=RK('fun',a,b,y0,h)
u3=0:0.001:h;
v3=Y2(1)+0*u3;
plot(u3,v3,'y--')
hold on
u4=0:0.001:5*h;
v4=Y2(5)+0*u4;
plot(u4,v4,'y--')
hold on
v5=0:0.001:Y2(1);
u5=h+0*v5;
plot(u5,v5,'k--')
hold on
v6=0:0.001:Y2(5);
u6=5*h+0*v6;
plot(u6,v6,'k--') hold on Y2=[d,Y2]; scatter(t,Y2,'r') hold on plot(t,Y2) title('¾«È·½âÇúÏßÓë½üËƽâÕÛÏß','fontsize',10,'fontweight','bold') text(0.735,0.7,'\leftarrowy=Euler·¨½üËƽâÕÛÏß','fontsize',8)
text(0.15,0.775,'\leftarrowy=¾-µä4½×R
unge-Kutta·¨½üËƽâÕÛÏß','fontsize',8)
x=[0.1,0.5];
y=exp(-2*x)+x.^2;
e1=abs(y(1)-Y1(2))
e2=abs(y(1)-Y2(2))
e1=abs(y(2)-Y1(6))
e2=abs(y(2)-Y2(6))
Euler法程序： function Y=Euler(f,a,b,y0,h) m=(b-a)/h;Y=zeros(1,m);x=a;d=y0; for n=1:m K=feval(f,x,y0); x=x+h; y0=y0+h*K; Y(n)=y0; end 定义函数： function z=fun(t,y) z=-2*y+2*t^2+2*t; 经典4阶Runge—Kutta法程序：
function Y=RK(f,a,b,y0,h)
m=(b-a)/h;Y=zeros(1,m);x=a;d=y0;
for n=1:m
K1=feval(f,x,y0);
x=x+0.5*h;y1=y0+0.5*h*K1;
K2=feval(f,x,y1);
y2=y0+0.5*h*K2;
K3=feval(f,x,y2);
x=x+0.5*h;y3=y0+h*K3;
K4=feval(f,x,y3);
y0=y0+h*(K1+2*K2+2*K3+K4)/6;
Y(n)=y0;
end
出师表
两汉：诸葛亮
先帝创业未半而中道崩殂，今天下三分，益州疲弊，此诚危急存亡之秋也。然侍卫之臣
不懈于内，忠志之士忘身于外者，盖追先帝之殊遇，欲报之于陛下也。诚宜开张圣听，以光
先帝遗德，恢弘志士之气，不宜妄自菲薄，引喻失义，以塞忠谏之路也。
宫中府中，俱为一体；陟罚臧否，不宜异同。若有作奸犯科及为忠善者，宜付有司论其
刑赏，以昭陛下平明之理；不宜偏私，使内外异法也。
侍中、侍郎郭攸之、费祎、董允等，此皆良实，志虑忠纯，是以先帝简拔以遗陛下：愚
以为宫中之事，事无大小，悉以咨之，然后施行，必能裨补阙漏，有所广益。
将军向宠，性行淑均，晓畅军事，试用于昔日，先帝称之曰“能”，是以众议举宠为督：愚
以为营中之事，悉以咨之，必能使行阵和睦，优劣得所。
亲贤臣，远小人，此先汉所以兴隆也；亲小人，远贤臣，此后汉所以倾颓也。先帝在时，
每与臣论此事，未尝不叹息痛恨于桓、灵也。侍中、尚书、长史、参军，此悉贞良死节之臣，
愿陛下亲之、信之，则汉室之隆，可计日而待也。
臣本布衣，躬耕于南阳，苟全性命于乱世，不求闻达于诸侯。先帝不以臣卑鄙，猥自枉
屈，三顾臣于草庐之中，咨臣以当世之事，由是感激，遂许先帝以驱驰。后值倾覆，受任于
败军之际，奉命于危难之间，尔来二十有一年矣。
先帝知臣谨慎，故临崩寄臣以大事也。受命以来，夙夜忧叹，恐托付不效，以伤先帝之
明；故五月渡泸，深入不毛。今南方已定，兵甲已足，当奖率三军，北定中原，庶竭驽钝，
攘除奸凶，兴复汉室，还于旧都。此臣所以报先帝而忠陛下之职分也。至于斟酌损益，进尽
忠言，则攸之、祎、允之任也。
愿陛下托臣以讨贼兴复之效，不效，则治臣之罪，以告先帝之灵。若无兴德之言，则责
攸之、祎、允等之慢，以彰其咎；陛下亦宜自谋，以咨诹善道，察纳雅言，深追先帝遗诏。
臣不胜受恩感激。
今当远离，临表涕零，不知所言。
