四阶龙格库塔和向前欧拉方法和隐式欧拉法以及改进欧拉法教材
6.1欧拉方法与龙格-库塔法
与相比较得 : ω1 + ω 2 = 1 1 αω 2 = 2 βω 2 = 1 2 由于四个参数，三个方 程，因此有一个自由参 数， 即解答不唯一。
1 1 (1) 取ω1 = , 可得ω2 = , α = β = 1, 此时算式为 2 2 1 yn +1 = yn + 2 (k1 + k 2 ) k1 = hf ( xn , yn ) k = hf ( x + h, y + k ) n n 1 2 这是改进的Euler方法。
常微分方程的解是一个函数，但是，计算机 没有办法对函数进行运算。因此，常微分方程的 数值解并不是求函数的近似，而是求解函数在某 些节点的近似值。 这种方法 ，称为数值离散方法。求的是在一 系列离散点列上，求未知函数y在这些点上的值的 近似。
6.1.1 欧拉(Euler)方法
设一阶常微分方程初值问题 dy = f ( x, y ) dx y ( x0 ) = y0 记 xn = a + nh (n = 0,1,2,...), h > 0为步长，一般总 假定为常数。该式的数值解是指通过某种方法 去获得解y ( x)在点xn 上的近似值yn , 即 y ( xn ) ≈ y n (n = 0,1,2,...)
= hf ( xn , yn ) + αh 2 f x′( xn , yn ) + β hk1 f y′ ( xn , yn ) + O(h 3 ) = hy′( xn ) + h 2 (αf x′( xn , yn ) + β f ( xn , yn ) f y′ ( xn , yn )) + O(h 3 )
于是 yn +1 = y ( xn ) + ω1hy′( xn ) + ω2 [hy′( xn ) + αh 2 f x′( xn , yn ) + βh 2 f ( xn , yn ) f y′ ( xn , yn ) + O(h 3 )] = y ( xn ) + (ω1 + ω2 )hy′( xn ) + [αω 2 f x′( xn , yn ) + βω 2 f ( xn , yn ) f y′ ( xn , yn )]h 2 + O(h 3 )
常微分方程数值解法-欧拉法、改进欧拉法与四阶龙格库塔法常微分方程数值解法
y( xn1)
y( xn
Байду номын сангаас
h)
y(xn )
hy'( xn )
h2 2!
y''( )
进一步: 令
h2 y( xn ) hy'( xn ) 2! y''( xn )
常微分方 yn1 y( xn1 ) , yn y( xn )
程数值解
法-欧拉法 yn1 yn hf ( xn , yn ) h2
、改进欧 y( xn1 ) yn1
2
max y''( x)
a xb
拉法和四
三、Euler方法
已 知 初 值 问 题 的 一 般 形式 为:
dy
dx
f (x, y)
a xb
(1)
y( x0 ) y0
常微分方 用差商近似导数 程数值解 问题转化为
yn1 yn dy
h
dx
法-欧拉法 yn1 yn hf ( xn , yn )
法-欧 y(拉0) 法1
、改进欧
拉法和四
四、几何意义
由 x0 , y0 出发取解曲线 y yx 的切线（存在！），则斜率
dy
f x0, y0
dx x y
,
0
0
常微分方 由于 f x0, y0 及 x0, y0 已知，必有切线方程。
由点斜式写出切程线方数程：值解
法、-改欧进拉欧法 ddxy y y0 x x0
常微分方 程数值解 能用解析方法求出精确解的微分方程为数不多,
而且有的方程即使有解析解,也可能由于解的表达
法-欧拉法 式非常复杂而不易计算,因此有必要研究微分方程
欧拉法与龙格库塔法比较分析
欧拉法与龙格库塔法比较分析解微分方程的欧拉法，龙格-库塔法简单实例比较欧拉方法(Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前EULER法、后退EULER法、改进的EULER法。
缺点:欧拉法简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。
因此欧拉格式一般不用于实际计算。
改进欧拉格式(向前欧拉公式):为提高精度，需要在欧拉格式的基础上进行改进。
采用区间两端的斜率的平均值作为直线方程的斜率。
改进欧拉法的精度为二阶。
算法:微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值。
对于常微分方程:dy,fxy(,)xab,[,] dxyay(), 0'可以将区间分成段，那么方程在第点有，再用n[,]abyxfxyx()(,()),xiiii 向前差商近似代替导数则为:((1)())yxyx，,ii,fxyx(,()) iihh在这里，是步长，即相邻两个结点间的距离。
因此可以根据点和的数xyii值计算出来: yi，1yyhfxy,，，(,)iL,0,1,2,？ iiii，1这就是向前欧拉公式。
改进的欧拉公式:将向前欧拉公式中的导数改为微元两端导数的平均，即上式便是梯形的fxy(,)ii欧拉公式。
可见，上式是隐式格式，需要迭代求解。
为了便于求解，使用改进的欧拉公式: 数值分析中，龙格,库塔法(Runge-Kutta)是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。
实际上，龙格-库塔法是欧拉方法的一种推广，向前欧拉公式将导数项简单取为,而改进的欧拉公式将导数项取为两端导数fxy(,)nn 的平均。
龙格-库塔方法的基本思想:在区间内多取几个点，将他们的斜率加权平均，作为导数的近似。
[,]xxnn，1龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。
令初值问题表述如下。
'yty(),yfty,(,) 00则，对于该问题的RK4由如下方程给出:h,，，，，(22)yykkkk nn，112346其中kfty,(,) 1nnhh(,)kftyk,，， 21nn22hh(,)kftyk,，， 32nn22kfthyhk,，，(,) 43nnh这样，下一个值由现在的值加上时间间隔和一个估算的斜率的乘积决yyn，1n定。
龙格库塔法介绍
yn
hf
(xn, yn ))]，
(x, y,h) 1[ f (x, y) f (x h, y hf (x, y))],
2
|
( x,
y1,
h)
(x,
y2 ,
h)
|
[L
2
L 2
(1
hL)]
|
y1
y2
|,
L
L(1
h0L)，h 2
h0.
类似地,不难验证其他龙格 库塔方法的收敛性.
这里c1,c2,c3,2,3, 21, 31, 32均为待定参数.
Tn1 y(xn1) yn1 O(h4 )
(3.11)
c1 c2 c3 1
2
21
3 31 32
c22
c33
1 2
cc232223c2332
将步长折半，从xn用两步求xn1处的近似值，则有
y(xn1)
h
yn21
2c
h 2
5
.
从而
h
y ( xn 1) y ( xn 1)
yn21 ynh1
1， 16
得到事后估计式：
y ( xn 1)
h
yn21
1 15
(
h
yn21
ynh1).
通过检查步长折半前后计算结果的偏差，
y(x) (x, y(x),0) 0 p 1 单步法(4.1)收敛. 定义4 若单步法(4.1)增量函数(x, y,h)是否满足
第四讲龙格-库塔方法
龙格-库塔方法3.2 Runge-Kutta法3.2.1 显式Runge-Kutta法的一般形式上节已给出与初值问题(1.2.1)等价的积分形式(3.2.1)只要对右端积分用不同的数值求积公式近似就可得到不同的求解初值问题(1.2.1)的数值方法，若用显式单步法(3.2.2)当，即数值求积用左矩形公式，它就是Euler法(3.1.2)，方法只有一阶精度，若取(3.2.3)就是改进Euler法，这时数值求积公式是梯形公式的一种近似，计算时要用二个右端函数f的值，但方法是二阶精度的.若要得到更高阶的公式，则求积分时必须用更多的f值，根据数值积分公式，可将(3.2.1)右端积分表示为注意，右端f中还不能直接得到，需要像改进Euler法(3.1.11)一样，用前面已算得的f值表示为(3.2.3)，一般情况可将(3.2.2)的 表示为(3.2.4)其中这里均为待定常数，公式(3.2.2)，(3.2.4)称为r级的显式Runge-Kutta法，简称R-K方法.它每步计算r个f值(即)，而k由前面(i-1)个已算出的表示，故公式是显式的.例i如当r=2时，公式可表示为(3.2.5) 其中.改进Euler 法(3.1.11)就是一个二级显式R-K 方法.参数取不同的值，可得到不同公式.3.2.2 二、三级显式R-K 方法对r=2的显式R-K 方法(3.2.5)，要求选择参数，使公式的精度阶p 尽量高，由局部截断误差定义11122211()()[(,())(,)]n n n n n n n T y x y x h c f x y x c f x a h y b hk ++=--+++ (3.2.6) 令，对(3.2.6)式在处按Taylor 公式展开，由于将上述结果代入(3.2.6)得要使公式(3.2.5)具有的阶p=2,即，必须(3.2.7)即由此三式求的解不唯一.因r=2，由（3.2.5）式可知，于是有解(3.2.8)它表明使(3.2.5)具有二阶的方法很多，只要都可得到二阶精度R-K方法.若取，则，则得改进Euler法(3.1.11)，若取，则得，此时(3.2.5)为(3.2.9)其中称为中点公式.改进的Euler法(3.1.11)及中点公式(3.2.9)是两个常用的二级R-K方法，注意二级R-K方法只能达到二阶，而不可能达到三阶.因为r=2只有4个参数，要达到p=3则在(3.2.6)的展开式中要增加3项，即增加三个方程,加上(3.2.7)的三个方程，共计六个方程求4个待定参数，验证得出是无解的.当然r=2,p=2的R-K方法(3.2.5)当取其他数时，也可得到其他公式，但系数较复杂，一般不再给出.对r=3的情形，要计算三个k值，即其中将按二元函数在处按Taylor公式展开，然后代入局部截断误差表达式，可得可得三阶方法，其系数共有8个，所应满足的方程为(3.2.10)这是8个未知数6个方程的方程组，解也是不唯一的，通常.一种常见的三级三阶R-K方法是下面的三级Kutta方法：(3.2.11)附：R-K 的三级Kutta 方法程序如下function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h;y = zeros(N+1,1); y(1) = y0; x = a:h:b;var = findsym(f); for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6; %满足c1+c2+c3=1,(1/6 4/6 1/6)endformat short; 3.2.3 四阶R-K 方法及步长的自动选择利用二元函数Taylor 展开式可以确定(3.2.4)中r=4,p=4的R-K 方法，其迭代公式为111223344()n n y y h c k c k c k c k +=++++其中1(,)n n k f x y =，2221(,(,))n n n n k f x a h y b hf x y =++，而33311322(,)n n k f x a h y b hk b hk =+++ 44411422433(,)n n k f x a h y b hk b hk b hk =++++共计13个参数待定，Taylor 展开分析局部截断误差，使得精度达到四阶，即误差为5()O h 。
四阶龙格——库塔法
四阶龙格——库塔法2013-2014（1）专业课程实践论文题目：四阶龙格—库塔法2'''()11()()()......()()2!!pp p n n n n n h h y y x hy x y x y x O h p ++=+++++ 进行离散化，所得计算公式必为p 阶方法，式中'''''()(,),()(,)(,)(,)....x y x f x y y x f x y f x y f x y ==++由此，我们能够想到，通过提高泰勒展开式的阶数，可以得到高精度的数值方法，从理论上讲，只要微分方程的解()y x 充分光滑，泰勒展开方法可以构造任意的有限阶的计算公式，但事实上，具体构造这种公式往往相当困难，因为符合函数(,())f x y x 的高阶导数常常是很烦琐的，因此，泰勒展开方法一般不直接使用，但是我们可以间接使用泰勒展开方法，求得高精度的计算方法。
首先，我们对欧拉公式和改进欧拉公式的形式作进一步的分析。
如果将欧拉公式和改进的欧拉公式改写成如下的形式：欧拉公式{111(,)n n n n y y hK K f x y +==+改进的欧拉公式11211()22n n y y h K K +=++， 1(,)n n K f x y =，21(,)n n K f x h y hK =++。
这两组公式都是用函数(,)f x y 在某些点上的值的线性组合来计算1()n y x +的近似值1n y +，欧拉公式每前进一步，就计算一次(,)f x y 的值。
另一方面它是1()n y x +在n x 处的一阶泰勒展开式，因而是一阶方法。
改进的欧拉公式每前进一步，需要计算两次(,)f x y 的值。
另一方面它在(,)n n x y 处的泰勒展开式与1()n y x +在n x 处的泰勒展开式的前三项完全相同，因而是二阶方法。
这启发我们考虑用函数(,)f x y 在若干点上的函数值的线性组合来构造计算公式。
常微分方程的数值解法(欧拉法、改进欧拉法、泰勒方法和龙格库塔法)
［例1］用欧拉方法与改进的欧拉方法求初值问题h 的数值解。
在区间[0,1]上取0.1［解］欧拉方法的计算公式为x0=0;y0=1;x(1)=0.1;y(1)=y0+0.1*2*x0/(3*y0^2);for n=1:9x(n+1)=0.1*(n+1);y(n+1)=y(n)+0.1*2*x(n)/(3*y(n)^2);end;xy结果为x =Columns 1 through 80.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 9 through 100.9000 1.0000y =Columns 1 through 81.0000 1.0067 1.0198 1.0391 1.0638 1.0932 1.1267 1.1634 Columns 9 through 101.2028 1.2443改进的欧拉方法其计算公式为本题的精确解为()y x=x0=0;y0=1;ya(1)=y0+0.1*2*x0/(3*y0^2);y(1)=y0+0.05*(2*x0/(3*y0^2)+2*x0/(3*ya^2));for n=1:9x(n+1)=0.1*(n+1);ya(n+1)=ya(n)+0.1*2*x(n)/(3*ya(n)^2);y(n+1)=y(n)+0.05*(2*x(n)/(3*y(n)^2)+2*x(n+1)/(3*ya(n+1)^2));end;xy结果为x =Columns 1 through 80.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 9 through 100.9000 1.0000y =Columns 1 through 81.0000 1.0099 1.0261 1.0479 1.0748 1.1059 1.1407 1.1783 Columns 9 through 101.2183 1.2600［例2］用泰勒方法解x＝0.1, 0.2, …, 1.0处的数值解，并与精确解进行比较。
用改进的欧拉方法和四阶龙格库塔方法解初值问题
用改进的欧拉方法和四阶龙格-库塔方法解初值问题一、题目：取步长2.0=h ，分别用改进的欧拉方法和四阶龙格—库塔方法解初值问题⎩⎨⎧=≤≤+=.1)0(;10,'y x y x y 并比较结果。
二、基本思想：1.改进的欧拉方法改进的欧拉方法用梯形公式计算()()()()+1+1-=,.n n x n n x y x y x f x y x dx ⎰中的积分，并以n y 和+1n y 分别表示()n y x 和()+1n y x 的近似值，就得到()()()+1+1+1=+,+,.2n n n n n n h y y f x y f x y （梯形公式） 梯形公式也是一个一步法公式。
由于公式右端也含有未知的+1n y ，故被称作是隐式的。
隐式格式实际上是关于+1n y 的一个函数方程。
为了避免解方程，可以采用欧拉方法计算初始值，再由梯形公式计算。
这样建立起来的计算格式称为改进的欧拉格式：()()()()+1+1+1=+,,=+,+,.2p n n n n n n n n n y y hf x y h y y f x y f x y ⎧⎪⎨⎪⎩ 梯形公式也可以采用迭代法求解。
如果仍然采用欧拉方法计算迭代初值，那么计算格式就是[]()[]()[]()()0+1+1+1+1+1=+,,=+,+,.2n nn n k k n n n n n n y y hf x y h y y f x y f x y ⎧⎪⎨⎪⎩ 由于已假定(),f x y 满足里普希兹条件，所以有[][][]()[]()[][]+1-1-1+1+1+1,+1+1+1+1+1-=-,-.22k k k k k k n n n n n n n n h hL y y f x y f x y y y ≤ 从而，迭代的收敛条件是0<<1.2hL 2.四阶龙格-库塔方法龙格—库塔方法不是用求导数的办法，而是用计算不同点上()y x f ,的函数值，然后对这些函数值作线性拟合，构造近似公式。
龙格-库塔方法基本原理3
x
2019/12/25
12
二阶龙格—库塔法
在[xi, xi+1]上取两点xi和xi+a2= xi +a2h,以该两点处的斜率值K1和K2的 加权平均(或称为线性组合)来求取平均斜率k*的近似值K，即
K c1 K 1 c 2 K 2
式中:K1为xi点处的切线斜率值 K1 =hf(xi, yi)=hy'(xi) K2为xi +a2h点处的切线斜率值,比照改进的欧拉法,将xi+a2视为xi+1，
R-K方法不是直接使用Taylor级数,而是利用它的思想
2019/12/25
5
Runge-Kutta 方法是一种高精度的单步法,简称R-K法
龙格-库塔(R-K)法的基本思想 Euler公式可改写成
yi 1 yi K K hf ( x i , y i )
则yi+1的表达式与y(xi+1)的Taylor展开式的前两项完全相同,即局部截 断误差为O(h2)。
即
y ( x i 1 ) y ( x i ) h y ( i )
2019/12/25
10
引入记号
y ( xi1 ) y ( xi ) K
K h y ( i ) hf i , y ( i )
yi 1 yi K
K 可以认为是
y y ( x ) 在区间 [ x i , x i 1 ]上的平均斜率
假设式
y' =f(x,y) (a≤x≤b)
中的 f(x,y) 充分光滑,将y(xi+1)在x i点作Taylor展开，若取右端不同的有限 项作为y(xi+1)的近似值，就可得到计算y(xi+1)的各种不同截断误差的数值 公式。
欧拉法、改进的欧拉法、龙格-库塔法求解初值问题
欧拉法、改进的欧拉法、龙格-库塔法求解初值问题欧拉法、改进的欧拉法、龙格-库塔法求解初值问题简介通过求解简单的初值问题：dudx =f (x ,u )(1)u (x 0)=u 0(2)引⼊欧拉法、改进的欧拉法、龙格-库塔法等。
前期准备数值解法的基本思想就是先对x 和u(x)在区间[x0,∞)上进⾏离散化，然后构造递推公式，再进⼀步得到u(x)u(x) u(x)u(x)u(x)u(x)在这些位置的近似取值。
取定步长h ，令x n =x 0+nh (n =±1,±2,⋯)得到离散的位置：x 1,x 2,⋯,x n ,u(x)在这些点精确取值为：u (x 1),u (x 2),⋯,u (x n )利⽤数值解法得到的这些点的近似取值，u 1,u 2,⋯,u n欧拉法欧拉法的核⼼就是将导数近似为差商。
将导数近似为向前差商，则有：du dxx =x n≈u x n +1−u x nh代⼊(1)式，有：u x n +1=y x n +hf x n ‖u x n⽤u n +1和 u n 代替u (x n +1)和u (x n )，得：u n +1=u n +hf x n ,u n因此，若知道u 0我们就可以递归出u 1,u 2,⋯如果将导数近似为向后差商:du dxx =x n≈u x n −u x n −1h类似的，就可以得到：u n −1=u n −hf x n ,u n这样，若知道u 0我们就可以递归出u −1,u −2⋯改进的欧拉法对(1)式在[x n ,x n +1]上积分，可得：u x n +1=u x n +∫xn +1x nf (x ,u )dx其中，n =0,1,⋯⽤不同⽅式来近似上式的积分运算，就会得到不同的递推公式。
若使⽤左端点计算矩形⾯积并取近似：∫x n +1x nf (x ,u )dx ≈hf x n +1,u x n +1代⼊上式得：{|()()()()(())()|()()()()()(())u n +1=u n +hf (x n ,u n )若使⽤梯形的⾯积做近似：∫x n +1x nf (x ,y )dx ≈h2f x n ,u x n+f x n +1,ux n +1得到：u n +1=u n +h2f x n ,u n +f x n +1,u n +1欧拉法虽然精度偏低，但它是显式的，可直接得到结果。
四阶龙格——库塔法
2013-2014（1）专业课程实践论文题目：四阶龙格—库塔法一、算法理论由定义可知，一种数值方法的精度与局部截断误差()po h 有关，用一阶泰勒展开式近似函数得到欧拉方法，其局部截断误差为一阶泰勒余项2()o h ，故是一阶方法，完全类似地若用p 阶泰勒展开式2'''()11()()()......()()2!!pp p n n n n n h h y y x hy x y x y x O h p ++=+++++ 进行离散化，所得计算公式必为p 阶方法，式中'''''()(,),()(,)(,)(,)....x y x f x y y x f x y f x y f x y ==++由此，我们能够想到，通过提高泰勒展开式的阶数，可以得到高精度的数值方法，从理论上讲，只要微分方程的解()y x 充分光滑，泰勒展开方法可以构造任意的有限阶的计算公式，但事实上，具体构造这种公式往往相当困难，因为符合函数(,())f x y x 的高阶导数常常是很烦琐的，因此，泰勒展开方法一般不直接使用，但是我们可以间接使用泰勒展开方法，求得高精度的计算方法。
首先，我们对欧拉公式和改进欧拉公式的形式作进一步的分析。
如果将欧拉公式和改进的欧拉公式改写成如下的形式：欧拉公式{111(,)n n n n y y hK K f x y +==+改进的欧拉公式11211()22n n y y h K K +=++， 1(,)n n K f x y =，21(,)n n K f x h y hK =++。
这两组公式都是用函数(,)f x y 在某些点上的值的线性组合来计算1()n y x +的近似值1n y +，欧拉公式每前进一步，就计算一次(,)f x y 的值。
另一方面它是1()n y x +在n x 处的一阶泰勒展开式，因而是一阶方法。
改进的欧拉公式每前进一步，需要计算两次(,)f x y 的值。
用Euler法、改进的Euler法及四阶的龙格库塔法求解初值问题
微分方程数值解课程设计题目1（30分）分别用Euler 法、改进的Euler 法及四阶的龙格库塔法求解初值问题：⎪⎩⎪⎨⎧=-=1)0(2'u ut u u 10≤<th 分别取0.1和0.2，要求计算并绘出图形，最后比较三种算法的精度。
（1）.先构建初值问题的函数M 文件： function z=fun(x,y)z=y-2*x/y; （2）.Euler 法 M 文件：function E=Euler(fun,x0,y0,h,N) x=zeros(1,N+1);y=zeros(1,N+1); x(1)=x0;y(1)=y0;for n=1:Nx(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(fun,x(n),y(n)); end T=[x;y]① .当h=0.1时>> Euler('fun',0,1,0.1,10) T = Columns 1 through 100 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.90001.0000 1.1000 1.1918 1.2774 1.3582 1.4351 1.5090 1.5803 1.6498 1.7178 Column 11 1.00001.7848②.当h=0.2时>> Euler('fun',0,1,0.2,5)T = 0 0.2000 0.4000 0.6000 0.8000 1.00001.0000 1.2000 1.3733 1.5315 1.6811 1.8269（3）.改进的Euler 法：M 文件：Euler_modify.mfunction E=Euler_modify(fun,x0,y0,h,N) x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=1:Nx(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n));y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n))+feval(fun,x(n+1),z0));endT=[x;y]①.当h=0.1时>> Euler_modify('fun',0,1,0.1,10)T = Columns 1 through 90 0.1000 0.2000 0.3000 0.4000 0.5000 0.60000.7000 0.80001.0000 1.0959 1.1841 1.2662 1.3434 1.4164 1.4860 1.5525 1.6165Columns 10 through 110.9000 1.00001.6782 1.7379②.当h=0.2时>> Euler_modify('fun',0,1,0.2,5)T = 0 0.2000 0.4000 0.6000 0.8000 1.00001.0000 1.1867 1.3483 1.4937 1.6279 1.7542（4）.四阶R_K方法：M文件：function [x,y]=Rk_N4(f,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=1:Nx(n+1)=x(n)+h;k1=h*feval(f,x(n),y(n));k2=h*feval(f,x(n)+1/2*h,y(n)+1/2*k1);k3=h*feval(f,x(n)+1/2*h,y(n)+1/2*k2);k4=h*feval(f,x(n)+h,y(n)+k3);y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);end①.当h=0.1时>> [x,y]=Rk_N4('fun',0,1,0.1,10)x = Columns 1 through 90 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000Columns 10 through 110.9000 1.0000y = Columns 1 through 91.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125Columns 10 through 111.6733 1.7321②.当h=0.2时>> [x,y]=Rk_N4('fun',0,1,0.2,5)x = 0 0.2000 0.4000 0.6000 0.8000 1.0000y = 1.0000 1.1832 1.3417 1.4833 1.6125 1.7321（5）作图：1): h=0.1>> x=0:0.1:1;>> y_euler=[1.0000 1.1000 1.1918 1.2774 1.3582 1.4351 1.5090 1.5803 1.6498 1.7178 1.7848];>> y_euler_modify=[1.0000 1.0959 1.1841 1.2662 1.3434 1.4164 1.4860 1.5525 1.6165 1.6782 1.7379];>> y_RK_N4=[1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321];plot(x,y_euler,'bo:',x,y_euler_modify,'r*-',x,y_RK_N4,'gv--');title('误差分析');xlabel('x轴');y label('y轴');text(0.3,1.3,'Euler法');text(0.4,1.35,'Euler改进法');text(0.5,1.4,'R-k法'); text(0.8,1.0 2,'作者：李靖');text(0.8,1.07,'日期：2010.6.19');text(0.4,1.75,'h=0.1');legend('Euler','Euler改进法','R_K法');grid on2): h=0.2>> x=0:0.2:1;>> y_euler=[1.0000 1.2000 1.3733 1.5315 1.6811 1.8269];>> y_euler_modify=[1.0000 1.1867 1.3483 1.4937 1.6279 1.7542];>> y_RK_N4=[ 1.0000 1.1832 1.3417 1.4833 1.6125 1.7321];plot(x,y_euler,'bo:',x,y_euler_modify,'r*-',x,y_RK_N4,'gv--');title('误差分析');xlabel('x轴');y label('y轴');text(0.3,1.3,'Euler法');text(0.4,1.35,'Euler改进法');text(0.5,1.4,'R-k法'); text(0.8,1.0 2,'作者：xx');text(0.8,1.07,'日期：2010.6.19'); text(0.4,1.75,'h=0.2');legend('Euler','Euler改进法', 'R_K法');grid on0.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8误差分析x 轴y 轴Euler 法Euler 改进法R-k 法作者：李靖日期：2010.6.19h=0.1EulerEuler 改进法R K 法00.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.81.9误差分析x 轴y 轴Euler 法Euler 改进法R-k 法作者：李靖日期：2010.6.19h=0.2EulerEuler 改进法R K 法2.（20分）编写一个程序用tylor 级数法求解问题：⎩⎨⎧==1)0('u ut u 10≤<t取tylor 级数法的截断误差为O(21h )，即要用u(t)，u ’(t)，…，u (20)(t)的值。
用欧拉法、改进欧拉法、四阶龙格-库塔法解初值问题
常微分方程的数值解一、实验名称：用欧拉法、改进欧拉法、四阶龙格-库塔法解初值问题[])( b a, x ),(00⎪⎩⎪⎨⎧=∈=y x y y x f dx dy 要求输出：x ，真值，欧拉法，改进欧拉法，龙格库塔法二、算法：（1）欧拉法公式：),(1n n n n y x f h y y +=+（2）改进欧拉法公式：[]),(),(2111+++++=n n n n n n y x f y x f h y y （3）龙格—库塔法公式： ()⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=⎪⎭⎫ ⎝⎛++=⎪⎭⎫ ⎝⎛++==++++=+342312143211,2,212,21),()22(6hk y h x f k k h y h x f k k h y h x f k y x f k k k k k h y y n n n n n n n n n n 三、源程序（见附件）1. 常微分方程见f.m2. 欧拉法、改进欧拉法、龙格库塔法及真值对比见cwffc.m四、数值例子以讲义例9.2为例：⎪⎩⎪⎨⎧=-=1)0(2y y x y dxdy 在区间[0，1.5]上，取h=0.11.在命令窗口输入: cwffc(@f,1,0.1,0,1.5)、回车输出：欧拉法、改进欧拉法、龙格库塔法的计算结果及真值五、分析总结由以上的实例可以看出这三种方法的计算精度：龙格库塔>改进欧拉法>欧拉法，由于欧拉法采取的实际是左矩形公式，随着计算次数的增多，其误差越来越大，改进欧拉法较欧拉法的精度要高，而龙格库塔法具有收敛、稳定、计算过程中可以改变步长等优点。
六、操作手册1.双击启动matlab；2.在命令窗口依次输入cwffc(@fx,y0,h,a,b)，回车（输入时给y0赋初始值，h赋步长，a赋左端点值，b赋右端点值，常微分方程的修改在f.m 中进行），即可输出真值和欧拉法、改进欧拉法、龙格库塔法的计算结果。
常微分方程的数值解法(欧拉法、改进欧拉法、泰勒方法和龙格库塔法)
［例1］用欧拉方法与改进的欧拉方法求初值问题h 的数值解。
在区间[0,1]上取0.1［解］欧拉方法的计算公式为x0=0;y0=1;x(1)=0.1;y(1)=y0+0.1*2*x0/(3*y0^2);for n=1:9x(n+1)=0.1*(n+1);y(n+1)=y(n)+0.1*2*x(n)/(3*y(n)^2);end;xy结果为x =Columns 1 through 80.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 9 through 100.9000 1.0000y =Columns 1 through 81.0000 1.0067 1.0198 1.0391 1.0638 1.0932 1.1267 1.1634 Columns 9 through 101.2028 1.2443改进的欧拉方法其计算公式为本题的精确解为()y x=x0=0;y0=1;ya(1)=y0+0.1*2*x0/(3*y0^2);y(1)=y0+0.05*(2*x0/(3*y0^2)+2*x0/(3*ya^2));for n=1:9x(n+1)=0.1*(n+1);ya(n+1)=ya(n)+0.1*2*x(n)/(3*ya(n)^2);y(n+1)=y(n)+0.05*(2*x(n)/(3*y(n)^2)+2*x(n+1)/(3*ya(n+1)^2));end;xy结果为x =Columns 1 through 80.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 9 through 100.9000 1.0000y =Columns 1 through 81.0000 1.0099 1.0261 1.0479 1.0748 1.1059 1.1407 1.1783 Columns 9 through 101.2183 1.2600［例2］用泰勒方法解x＝0.1, 0.2, …, 1.0处的数值解，并与精确解进行比较。
龙格库塔法和欧拉法求解微分方程的比较
龙格库塔法和欧拉法求解微分方程的比较龙格库塔法和欧拉法是求解微分方程常用的数值方法。
它们都是通过离散化来近似求解连续的微分方程，但在精度和稳定性等方面存在一些差异。
我们来看一下欧拉法。
欧拉法是一种简单的数值方法，它通过将微分方程中的导数进行离散化来近似求解。
具体而言，欧拉法使用泰勒展开式将导数近似为差商，然后通过迭代的方式求解微分方程。
欧拉法的优点是简单易懂，计算速度快，但它的缺点也很明显，即精度低。
因为欧拉法只使用了一阶导数的信息，忽略了高阶导数的贡献，所以它的近似误差较大。
相比之下，龙格库塔法是一种更为精确的数值方法。
它通过对微分方程的积分进行逐步逼近，从而得到更精确的解。
龙格库塔法的基本思想是利用不同阶的差商来近似积分，从而得到更高阶的逼近解。
龙格库塔法的优点是精度高，收敛速度快，适用于各种类型的微分方程。
但是，龙格库塔法的计算量相对较大，需要进行多次迭代计算，所以在实际应用中可能会耗费一定的时间和计算资源。
总的来说，龙格库塔法相对于欧拉法来说是一种更为精确和稳定的数值方法。
它能够提供更精确的解，并且适用于各种类型的微分方程。
然而，在一些简单的问题中，欧拉法可能仍然是一个可行的选择，因为它简单易懂，计算速度快。
当我们在实际问题中遇到微分方程时，我们可以根据问题的特点选择合适的数值方法来求解。
如果我们追求更高的精度和稳定性，那么龙格库塔法是一个不错的选择。
但如果我们对精度要求不高，或者问题相对简单，那么欧拉法也可以作为一个快速且简单的解决方案。
需要指出的是，除了龙格库塔法和欧拉法，还有许多其他的数值方法可以用于求解微分方程。
每种方法都有其优缺点，我们需要根据具体的问题和求解要求来选择合适的方法。
在实际应用中，我们也可以结合多种方法来求解微分方程，以提高求解的精度和效率。
龙格库塔法和欧拉法是求解微分方程常用的数值方法。
它们在精度和稳定性等方面存在一些差异，我们可以根据具体问题的特点选择合适的方法来求解。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验九欧拉方法信息与计算科学金融崔振威201002034031一、实验目的：1、掌握欧拉算法设计及程序实现二、实验内容：1、p364-9.2.4、p386-9.5.6三、实验要求：主程序：欧拉方法（前项）：function [x,y]=DEEuler(f,a,b,y0,n);%f:一阶常微分方程的一般表达式的右端函数%a:自变量的取值下限%b:自变量的取值上限%y0：函数的初值%n:积分的步数if nargin<5,n=50;endh=(b-a)/n;x(1)=a;y(1)=y0;for i=1:nx(i+1)=x(i)+h;y(i+1)=y(i)+h*feval(f,x(i),y(i));endformat short欧拉方法（后项）：function [x,y]=BAEuler(f,a,b,y0,n);%f:一阶常微分方程的一般表达式的右端函数%a:自变量的取值下限%b:自变量的取值上限%y0：函数的初值%n:积分的步数if nargin<5,n=50;endh=(b-a)/n;x(1)=a;y(1)=y0;for i=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(f,x(i),y(i));y(i+1)=y(i)+h*feval(f,x(i+1),y1);endformat short梯形算法：function [I,step,h2] = CombineTraprl(f,a,b,eps)%f 被积函数%a,b 积分上下限%eps 精度%I 积分结果%step 积分的子区间数if(nargin ==3)eps=1.0e-4;endn=1;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;while abs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;for i=0:n-1x=a+h*i;x1=x+h;I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1));endendI=I2;step=n;h2=(b-a)/n;改进欧拉方法：function [x,y]=MoEuler(f,a,b,y0,n);%f:一阶常微分方程的一般表达式的右端函数%a:自变量的取值下限%b:自变量的取值上限%y0：函数的初值%n:积分的步数if nargin<5,n=50;endh=(b-a)/n;x(1)=a;y(1)=y0; for i=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(f,x(i),y(i)); y2=y(i)+h*feval(f,x(i+1),y1); y(i+1)=(y1+y2)/2; endformat short四阶龙格-库塔法：function y = DELGKT4_lungkuta(f, h,a,b,y0,varvec) %f:一阶常微分方程的一般表达式的右端函数 %h:积分步长%a ：自变量的取值下限 %b:自变量的取值上限%varvec ：常微分方程的变量组 format long; N = (b-a)/h;y = zeros(N+1,1); y(1) = y0; x = a:h:b;var = findsym(f); for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);%Funval 为程序所需要的函数 K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K2*h/2]); K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*K3]); y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6; endformat short;p364-9.2.4欧拉方法（前项）：1、22)(,1)0(,22'+-+-==-=-t t e t y y y t y t解:执行20步时：编写函数文件doty.m 如下： function f=doty(x,y); f=x^2-y;在Matlab 命令窗口输入：>> [x1,y1]=DEEuler('doty',0,2,1,20)可得到结果：x1 =0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.70000.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.50001.6000 1.7000 1.8000 1.90002.0000y1 =1.0000 0.9000 0.8110 0.7339 0.6695 0.6186 0.5817 0.55950.5526 0.5613 0.5862 0.6276 0.6858 0.7612 0.8541 0.96471.0932 1.2399 1.4049 1.5884 1.7906在Matlab 命令窗口输入：>> y3=-exp(-x1)+x1.^2-2*x1+2求得解析解：y3 =1.0000 0.9052 0.8213 0.7492 0.6897 0.6435 0.6112 0.59340.5907 0.6034 0.6321 0.6771 0.7388 0.8175 0.9134 1.02691.1581 1.3073 1.4747 1.6604 1.8647输入：>> plot(x1,y1,'o');hold on>> plot(x1,y3,'*');hold on可得近似值与解析解的图像比较：执行40步时：在Matlab 命令窗口输入：>> [x1,y1]=DEEuler('doty',0,2,1,40)可得到结果：x1 =0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.75000.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.15001.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.55001.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.95002.0000y1 =1.0000 0.9500 0.9026 0.8580 0.8162 0.7774 0.7417 0.7091 0.6798 0.6538 0.6312 0.6121 0.5967 0.5848 0.5767 0.5724 0.5719 0.5753 0.5826 0.5940 0.6094 0.6290 0.6526 0.68050.7126 0.7490 0.7897 0.8347 0.8841 0.9379 0.9961 1.05881.1260 1.1977 1.2739 1.3547 1.4401 1.5301 1.6247 1.7240 1.8279在Matlab 命令窗口输入：>> y3=-exp(-x1)+x1.^2-2*x1+2求得解析解：y3 =1.0000 0.9513 0.9052 0.8618 0.8213 0.7837 0.7492 0.7178 0.6897 0.6649 0.6435 0.6256 0.6112 0.6005 0.5934 0.5901 0.5907 0.5951 0.6034 0.6158 0.6321 0.6526 0.6771 0.70590.7388 0.7760 0.8175 0.8633 0.9134 0.9679 1.0269 1.09031.1581 1.2305 1.3073 1.3887 1.4747 1.5653 1.6604 1.7602 1.8647输入：>> plot(x1,y1,'o');hold on>> plot(x1,y3,'*');hold on可得近似值与解析解的图像比较：从上面结果可以看出，执行40步时近似解的值要接近于解析解，误差更小，结果更精确。
5、)1/(1)(,1)0(,222't t y y ty y -===解:执行20步时：编写函数文件doty.m 如下： function f=doty(x,y); f=2*x*y^2;在Matlab 命令窗口输入：>> [x1,y1]=DEEuler('doty',0,2,1,20) 可得到结果： x1 =0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 y1 =1.0e+176 *0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 6.4573在Matlab 命令窗口输入： >> y3=1./(1-x1.^2)求得解析解：y3 =1.0e+015 *0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.5036 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000输入：>> plot(x1,y1,'o');hold on>> plot(x1,y3,'*');hold on可得近似值与解析解的图像比较：执行40步时：在Matlab 命令窗口输入：>> [x1,y1]=DEEuler('doty',0,2,1,40)可得到结果：x1 =0.0000 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.75000.8000 0.8500 0.9000 0.9500 1.0000 1.05001.1000 1.1500 1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500 1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.95002.0000y1 =1.0e+224 *0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.00000.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3.6893 Inf Inf Inf Inf Inf Inf Inf Inf Inf 在Matlab 命令窗口输入： >> y3=1./(1-x1.^2) 求得解析解： y3 =1.0e+015 *0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -2.2518 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 输入：>> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on可得近似值与解析解的图像比较：欧拉方法（后项）：1、22)(,1)0(,22'+-+-==-=-t t e t y y y t y t解:执行20步时：编写函数文件doty.m 如下：function f=doty(x,y);f=x^2-y;在Matlab 命令窗口输入：>> [x1,y1]=BAEuler('doty',0,2,1,20)可得到结果：x1 =0.0000 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.70000.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.50001.6000 1.7000 1.8000 1.90002.0000y1 =1.0000 0.9110 0.8329 0.7665 0.7127 0.6719 0.6449 0.63230.6345 0.6520 0.6852 0.7345 0.8003 0.8829 0.9825 1.09951.2341 1.3864 1.5567 1.7452 1.9520在Matlab 命令窗口输入：>> y3=-exp(-x1)+x1.^2-2*x1+2求得解析解：y3 =1.0000 0.9052 0.8213 0.7492 0.6897 0.6435 0.6112 0.59340.5907 0.6034 0.6321 0.6771 0.7388 0.8175 0.9134 1.02691.1581 1.3073 1.4747 1.6604 1.8647输入：>> plot(x1,y1,'o');hold on>> plot(x1,y3,'*');hold on可得近似值与解析解的图像比较：执行40步时：在Matlab 命令窗口输入：>> [x1,y1]=BAEuler('doty',0,2,1,40)可得到结果：x1 =0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.75000.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.15001.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.55001.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.95002.0000y1 =1.0000 0.9526 0.9079 0.8658 0.8267 0.7904 0.7572 0.7272 0.7003 0.6768 0.6566 0.6399 0.6268 0.6172 0.6114 0.6092 0.6109 0.6164 0.6258 0.6392 0.6566 0.6780 0.7035 0.73320.7671 0.8052 0.8475 0.8942 0.9451 1.0005 1.0602 1.12431.1929 1.2660 1.3435 1.4256 1.5122 1.6034 1.6992 1.7996 1.9046在Matlab 命令窗口输入：>> y3=-exp(-x1)+x1.^2-2*x1+2求得解析解：y3 =1.0000 0.9513 0.9052 0.8618 0.8213 0.7837 0.7492 0.71780.6897 0.6649 0.6435 0.6256 0.6112 0.6005 0.5934 0.5901 0.5907 0.5951 0.6034 0.6158 0.6321 0.6526 0.6771 0.7059 0.7388 0.7760 0.8175 0.8633 0.9134 0.9679 1.0269 1.0903 1.1581 1.2305 1.3073 1.3887 1.4747 1.5653 1.6604 1.7602 1.864 输入：>> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on可得近似值与解析解的图像比较：从上面结果可以看出，执行40步时近似解的值要接近于解析解，误差更小，结果更精确。