东南大学_数值分析_第六章_常微分方程数值解法
- 格式:doc
- 大小:152.00 KB
- 文档页数:6
许多实际问题的数学模型是微分方程或微分方程的定解问题。
如物体运动、电路振荡、化学反映及生物群体的变化等。
常微分方程可分为线性、非线性、高阶方程与方程组等类;线性方程包含于非线性类中,高阶方程可化为一阶方程组。
若方程组中的所有未知量视作一个向量,则方程组可写成向量形式的单个方程。
因此研究一阶微分方程的初值问题⎪⎩⎪⎨⎧=≤≤=0)(),(y a y bx a y x f dxdy, (9-1) 的数值解法具有典型性。
常微分方程的解能用初等函数、特殊函数或它们的级数与积分表达的很少。
用解析方法只能求出线性常系数等特殊类型的方程的解。
对非线性方程来说,解析方法一般是无能为力的,即使某些解具有解析表达式,这个表达式也可能非常复杂而不便计算。
因此研究微分方程的数值解法是非常必要的。
只有保证问题(9-1)的解存在唯一的前提下,研究其数值解法或者说寻求其数值解才有意义。
由常微分方程的理论知,如果(9-1)中的),(y x f 满足条件(1)),(y x f 在区域} ),({+∞<<∞-≤≤=y b x a y x D ,上连续; (2)),(y x f 在上关于满足Lipschitz 条件,即存在常数,使得y y L y x f y x f -≤-),(),(则初值问题(9-1)在区间],[b a 上存在惟一的连续解)(x y y =。
在下面的讨论中,我们总假定方程满足以上两个条件。
所谓数值解法,就是求问题(9-1)的解)(x y y =在若干点b x x x x a N =<<<<= 210处的近似值),,2,1(N n y n =的方法。
),,2,1(N n y n =称为问题(9-1)的数值解,n n x x h -=+1称为由到1+n x 的步长。
今后如无特别说明,我们总假定步长为常量。
建立数值解法,首先要将微分方程离散化,一般采用以下几种方法: (1) 用差商近似导数在问题(9-1)中,若用向前差商hx y x y n n )()(1-+代替)(n x y ',则得)1,,1,0( ))(,()()(1-=≈-+N n x y x f hx y x y n n n n n)(n x y 用其近似值代替,所得结果作为)(1+n x y 的近似值,记为1+n y ,则有 1(,) (0,1,,1)n n n n y y hf x y n N +=+=-这样,问题(9-1)的近似解可通过求解下述问题100(,) (0,1,,1)()n n n n y y hf x y n N y y x +=+=-⎧⎨=⎩(9-2)得到,按式(9-2)由初值经过步迭代,可逐次算出N y y y ,,21。
《计算机数学基础》数值部分第五单元辅导14 常微分方程的数值解法一、重点内容 1. 欧拉公式:),...,,,(),()(1-210=⎩⎨⎧+=+=≈01+1+n k kh x x y x hf y y x y kk k k k k局部截断误差是O (h 2)。
2. 改进欧拉公式:预报-校正公式:⎪⎩⎪⎨⎧++=+=++++)],(),([2),(1111k k k k k k k k k k y x f y x f hy y y x hf y y 校正值预报值即 ))],(,(),([211k k k k k k k k y x hf y x f y x f hy y +++=++ 或表成平均的形式:⎪⎪⎪⎩⎪⎪⎪⎨⎧+21=+=+=1+1+)(),(),(c p k p k k c k k k p y y y y x hf y y y x hf y y改进欧拉法的局部截断误差是O (h 3)3. 龙格-库塔法二阶龙格-库塔法的局部截断误差是O (h 3) 三阶龙格-库塔法的局部截断误差是O (h 4) 四阶龙格−库塔法公式: )22(643211κκκκ++++=+hy y k k其中 κ1=f (x k ,y k );κ2=f (x n +12h ,y k +21h κ1);κ3=f (x k +12h ,y n +21h κ2);κ4=f (x k +h ,y k +h κ3)四阶龙格-库塔法的局部截断误差是O (h 5)。
二、实例例1 用欧拉法解初值问题⎩⎨⎧1=060≤≤0--='2)().(y x xy y y ,取步长h =0.2。
计算过程保留4位小数。
解h =0.2, f (x )=-y -xy 2。
首先建立欧拉迭代格式),,)((.),(210=-420=--=+=21+k y x y y hx hy y y x hf y y k k k kk k k k k k k当k =0,x 1=0.2时,已知x 0=0,y 0=1,有y (0.2)≈y 1=0.2×1(4-0×1)=0.8000当k =1,x 2=0.4时,已知x 1=0.2, y 1=0.8,有 y (0.4)≈y 2=0.2×0.8×(4-0.2×0.8)=0.614 4 当k =2,x 3=0.6时,已知x 2=0.4,y 2=0.6144,有 y (0.6)≈y 3=0.2×0.6144×(4-0.4×0.4613)=0.8000例2 用欧拉预报-校正公式求解初值问题⎩⎨⎧1=10=++'2)(sin y x y y y ,取步长h =0.2,计算y (0.2),y (0.4)的近似值,计算过程保留5位小数。
i.常微分方程初值问题数值解法i.1 常微分方程差分法考虑常微分方程初值问题:求函数()u t 满足(,), 0du f t u t T dt=<≤ (i.1a ) 0(0)u u = (i.1b)其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数,0u 和T 是给定的常数。
我们假设(,)f t u 对u 满足Lipschitz 条件,即存在常数L 使得121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-∀∈∈-∞∞ (i.2) 这一条件保证了(i.1)的解是适定的,即存在,唯一,而且连续依赖于初值0u 。
通常情况下,(i.1)的精确解不可能用简单的解析表达式给出,只能求近似解。
本章讨论常微分方程最常用的近似数值解法--差分方法。
先来讨论最简单的Euler 法。
为此,首先将求解区域[0,]T 离散化为若干个离散点:0110N N t t t t T -=<<<<= (i.3) 其中n t hn =,0h >称为步长。
在微积分课程中我们熟知,微商(即导数)是差商的极限。
反过来,差商就是微商的近似。
在0t t =处,在(i.1a )中用向前差商10()()u t u t h -代替微商du dt ,便得 10000()()(,())u t u t hf t u t ε=++如果忽略误差项0ε,再换个记号,用i u 代替()i u t 便得到1000(,)u u hf t u -=一般地,我们有1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-方法: (i.4) 从(i.1b) 给出的初始值0u 出发,由上式可以依次算出1,,N t t 上的差分解1,,N u u 。
下面我们用数值积分法重新导出 Euler 法以及其它几种方法。
为此,在区间1[,]n n t t +上积分常微分方程(i.1a ),得11()()(,())n n t n n t u t u t f t u t dt ++=+⎰ (i.5)用各种数值积分公式计算(i.5)中的积分,便导致各种不同的差分法。
常微分方程组数值解法一、引言常微分方程组是数学中的一个重要分支,它在物理、工程、生物等领域都有广泛应用。
对于一些复杂的常微分方程组,往往难以通过解析方法求解,这时候数值解法就显得尤为重要。
本文将介绍常微分方程组数值解法的相关内容。
二、数值解法的基本思想1.欧拉法欧拉法是最基础的数值解法之一,它的思想是将时间连续化,将微分方程转化为差分方程。
对于一个一阶常微分方程y'=f(x,y),其欧拉公式为:y_{n+1}=y_n+hf(x_n,y_n)其中h为步长,x_n和y_n为第n个时间点上x和y的取值。
2.改进欧拉法改进欧拉法是对欧拉法的改良,其公式如下:y_{n+1}=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},y_n+hf(x_n,y_n))] 3.四阶龙格-库塔方法四阶龙格-库塔方法是目前最常用的数值解法之一。
其公式如下:k_1=f(x_n,y_n)k_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}k_1)k_3=f(x_n+\frac{h}{2},y_n+\frac{h}{2}k_2)k_4=f(x_n+h,y_n+hk_3)y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)其中,k_i为中间变量。
三、常微分方程组的数值解法1.欧拉法对于一个二阶常微分方程组:\begin{cases} y'_1=f_1(x,y_1,y_2) \\ y'_2=f_2(x,y_1,y_2)\end{cases}其欧拉公式为:\begin{cases} y_{n+1,1}=y_{n,1}+hf_1(x_n,y_{n,1},y_{n,2}) \\y_{n+1,2}=y_{n,2}+hf_2(x_n,y_{n,1},y_{n,2}) \end{cases}其中,x_n和y_{n,i}(i=1, 2)为第n个时间点上x和y_i的取值。
第六章 常微分方程数值解法——RK 4法、AB 4法******(学号) *****(姓名)上机题目要求见教材P307,23题。
一、算法原理题目要求采用RK 4法和AB 4法求解最简单的常微分方程初值问题(,),()y f x y a x by a η'=≤≤⎧⎨=⎩ (1)为求解式(1),采用离散化方法,就是寻求解)(x y 在区间],[b a 上的一系列点<<<<<n x x x x 321上的近似值 ,,,,21n y y y 。
记1(1,2,)i i i h x x i -=-=表示相邻两个节点的间距,称为步长。
求微分方程数值解的主要问题:(1) 如何将微分方程(,)y f x y '=离散化,并建立求其数值解的递推公式; (2) 递推公式的局部截断误差、数值数n y 与精确解)(n x y 的误差估计; (3) 递推公式的稳定性与收敛性. a) Runge-Kutta 方法基本思想:通过在1[,]i i x x +多预报几个点求斜率,并将其加权平均作为k *的近似值,以此构造更高精度的计算公式。
如果每步计算四次函数 的值,完全类似的,可以导出局部截断误差为)(5h O 的四阶Runge-Kutta 公式(RK 4):1123412132431(22),6(,),(,),221(,),22(,).n n n n n n n n n n y y k k k k k f x y h h k f x y k h k f x h y k k f x h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ (2)b) Adams 显式公式Runge-Kutta 方法是单步法,计算1+n y 时,只用到n y , 而已知信息1-n y 、2-n y 等没有被直接利用。
可以设想如果充分利用已知信息1-n y ,2-n y ,…来计算1+n y ,那么不但有可能提高精度,而且大大减少了计算量,这就是构造所谓线性多步法的基本思想。
常微分方程的数值解法常微分方程是研究变量的变化率与其当前状态之间的关系的数学分支。
它在物理、工程、经济等领域有着广泛的应用。
解常微分方程的精确解往往十分困难甚至不可得,因此数值解法在实际问题中起到了重要的作用。
本文将介绍常见的常微分方程的数值解法,并比较其优缺点。
1. 欧拉方法欧拉方法是最简单的数值解法之一。
它基于近似替代的思想,将微分方程中的导数用差商近似表示。
具体步骤如下:(1)确定初始条件,即问题的初值。
(2)选择相应的步长h。
(3)根据微分方程的定义使用近似来计算下一个点的值。
欧拉方法的计算简单,但是由于误差累积,精度较低。
2. 改进欧拉方法为了提高欧拉方法的精度,改进欧拉方法应运而生。
改进欧拉方法通过使用两个点的斜率的平均值来计算下一个点的值。
具体步骤如下:(1)确定初始条件,即问题的初值。
(2)选择相应的步长h。
(3)根据微分方程的定义使用近似来计算下一个点的值。
改进欧拉方法相较于欧拉方法而言,精度更高。
3. 龙格-库塔法龙格-库塔法(Runge-Kutta)是常微分方程数值解法中最常用的方法之一。
它通过迭代逼近精确解,并在每一步中计算出多个斜率的加权平均值。
具体步骤如下:(1)确定初始条件,即问题的初值。
(2)选择相应的步长h。
(3)计算各阶导数的导数值。
(4)根据权重系数计算下一个点的值。
与欧拉方法和改进欧拉方法相比,龙格-库塔法的精度更高,但计算量也更大。
4. 亚当斯法亚当斯法(Adams)是一种多步法,它利用之前的解来近似下一个点的值。
具体步骤如下:(1)确定初始条件,即问题的初值。
(2)选择相应的步长h。
(3)通过隐式或显式的方式计算下一个点的值。
亚当斯法可以提高精度,并且比龙格-库塔法更加高效。
5. 多步法和多级法除了亚当斯法,还有其他的多步法和多级法可以用于解常微分方程。
多步法通过利用多个点的值来逼近解,从而提高精度。
而多级法则将步长进行分割,分别计算每个子问题的解,再进行组合得到整体解。
第六章常微分方程的数值解法第六章常微分方程的数值解法在自然科学研究和工程技术领域中,常常会遇到常微分方程的求解问题。
传统的数学分析方法仅能给出一些简单的、常系数的、经典的线性方程的解析表达式,不能处理复杂的、变系数的、非线性方程,对于这些方面的问题,只能求诸于近似解法和数值解法。
而且在许多实际问题中,确确实实并不总是需要精确的解析解,往往只需获得近似的解或者解在若干个点上的数值即可。
在高等数学课程中介绍过的级数解法和逐步逼近法,能够给出解的近似表达式,这一类方法称为近似解法。
还有一类方法是通过计算机来求解微分方程的数值解,给出解在一些离散点上的近似值,这一类方法称作为数值方法。
本章主要介绍常微分方程初值问题的数值解法,包括Euler 方法、Runge-Kutta 方法、线性多步法以及微分方程组与高阶微分方程的数值解法。
同时,对于求解常微分方程的边值问题中比较常用的打靶法与有限差分法作了一个简单的介绍。
§1 基本概念1.1 常微分方程初值问题的一般提法常微分方程初值问题的一般提法是求解满足如下条件的函数,,b x a x y ≤≤)(=<<=α)(),(a y bx a y x f dxdy, (1.1) 其中),(y x f 是已知函数,α是给定的数值。
通常假定上面所给出的函数),(y x f 在给定的区域},),{(+∞<≤≤=yb x a y x D 上面满足如下条件:(1) 函数),(y x f 在区域D 上面连续;(2) 函数),(y x f 在区域D 上关于变量y 满足Lipschitz(李普希茨)条件:212121,),(),(y y b x a y y L y x f y x f ?≤≤?≤?,, (1.2)其中常数L 称为Lipschitz(李普希茨)常数。
由常微分方程的基本理论可以知道,假如(1.1)中的),(y x f 满足上面两个条件,则常微分方程初值问题(1.1)对于任意给定的初始值α都存在着唯一的解,,b x a x y ≤≤)(并且该唯一解在区间[a,b]上是连续可微的。
第六章 常微分方程数值解法
——RK 4法、AB 4法
******(学号) *****(姓名)
上机题目要求见教材P307,23题。
一、算法原理
题目要求采用RK 4法和AB 4法求解最简单的常微分方程初值问题
(,),()y f x y a x b
y a η
'=≤≤⎧⎨
=⎩ (1)
为求解式(1),采用离散化方法,就是寻求解)(x y 在区间],[b a 上的一系列点
<<<<<n x x x x 321
上的近似值 ,,,,21n y y y 。
记1(1,2,)i i i h x x i -=-=表示相邻两个节点的间距,称为步长。
求微分方程数值解的主要问题:
(1) 如何将微分方程(,)y f x y '=离散化,并建立求其数值解的递推公式; (2) 递推公式的局部截断误差、数值数n y 与精确解)(n x y 的误差估计; (3) 递推公式的稳定性与收敛性. a) Runge-Kutta 方法
基本思想:通过在1[,]i i x x +多预报几个点求斜率,并将其加权平均作为k *的近似值,以此构造更高精度的计算公式。
如果每步计算四次函数 的值,完全类似的,可以导出局部截断误差为)(5h O 的四阶Runge-Kutta 公式(RK 4):
1
123412132431(22),6
(,),(,),
221(,),22(,).
n n n n n n n n n n y y k k k k k f x y h h k f x y k h k f x h y k k f x h y hk +⎧
=++++⎪⎪
=⎪⎪⎪
=++⎨⎪
⎪
=++⎪⎪=++⎪⎩ (2)
b) Adams 显式公式
Runge-Kutta 方法是单步法,计算1+n y 时,只用到n y , 而已知信息1-n y 、2
-n y 等没有被直接利用。
可以设想如果充分利用已知信息1-n y ,2-n y ,…来计算1+n y ,那么不但有可能提高精度,而且大大减少了计算量,这就是构造所谓线性多步法的基本思想。
我们知道常微分初值问题(1)与积分
1
1()()(,())d n n
x n n x y x y x f t y t t ++=+⎰
(3)
是等价的。
下面我们讨论最简单的多步法——Admas 方法。
为了导出求微分方程(1)的数值解法,我们将(3)右边积分的被积函数用插值多项式来近似替代,可以得到4阶精度的显式公式
11122331
[55(,)59(,)37(,)9(,)]
24
n n n n n n n n n n y y h f x y f x y f x y f x y +------=+-+- (4) 称为AB 4法。
二、流程图
1. 四阶Runge-Kutta 公式(RK 4)通用程序流程图1所示。
图1关于RK 4算法流程图
2.AB4法通用程序流程图2所示。
图2关于AB4算法流程图
三、计算代码
1.RK4完整代码
function y=Runge_Kutta4(f,y0,a,b,h)
%求解常微分方程初值问题
%Input - f is the function
% - a and b 定义域端点
% - y0 初始值
% - h 步长
%Output - y 输出
clc
y=y0;
ii=1;
for x=a+h:h:b
k1=feval(f,x,y(ii));
k2=feval(f,x+0.5*h,y(ii)+0.5*h*k1);
k3=feval(f,x+0.5*h,y(ii)+0.5*h*k2);
k4=feval(f,x+h,y(ii)+h*k3);
temp=y(ii)+h*(k1+2*k2+2*k3+k4)/6; y=[y,temp]; ii=ii+1; end
2. AB 4完整代码
function y=AB4(f,y0,a,b,h) %求解常微分方程初值问题 %Input - f is the function % - a and b 定义域端点 % - y0 初始值 % - h 步长 %Output - y 输出
y=Runge_Kutta4(f,y0,a,a+3*h,h); ii=4;
for x=a+3*h:h:b-h
temp=y(ii)+h*(55*feval(f,x,y(ii))-59*feval(f,x-h,y(ii-1))+... 37*feval(f,x-2*h,y(ii-2))-9*feval(f,x-3*h,y(ii-3)))/24; y=[y,temp]; ii=ii+1; end
四、计算结果及分析
对于初值问题
22 (0 1.5)
(0)3y x y x y '⎧=-≤≤⎨
=⎩
(5)
取步长0.1h =,采用RK4法和AB4法计算,并将计算结果与精确值
3()3/(1)y x x =+作比较,计算调用程序Runge_Kutta4和AB4,如下
%求解常微分方程dy/dx=-x^2*y^2的初值问题 %分别采用RK4法和AB4法求解,并与精确解作比较 clc,format long
f=inline('-x*x*y*y');%函数f(x,y)
y_rk4=Runge_Kutta4(f,3,0,1.5,0.1);%RK4法 disp('RK4法计算结果为:') fprintf('%8.5f ', y_rk4) fprintf('\n')
y_ab4=AB4(f,3,0,1.5,0.1);%AB4法 disp('AB4法计算结果为:') fprintf('%8.5f ', y_ab4) fprintf('\n') x=0:0.1:1.5;
y=3./(1+x.^3);%精确解
disp('精确结果为:')
fprintf('%8.5f ', y)
fprintf('\n')
error_rk4=abs(y-y_rk4);%RK4法误差
error_ab4=abs(y-y_ab4);%AB4法误差
disp('RK4法的误差为:')
fprintf('%8.5f ', error_rk4)
fprintf('\n')
disp('AB4法的误差为:')
fprintf('%8.5f ', error_ab4)
fprintf('\n')
在MATLAB中运行以上程序,计算结果如表1所示
表1 RK4法与AB4法计算结果及其误差
从表1可以看出,Rung-Kula方法和AB4法均能很好的实现常微分方程初值问题的解法,且具有较高的精度和较低的运算复杂度。
同时,根据计算的误差可见,RK4法较之于AB4法具有更高的计算精度,但计算量比AB4法稍微多一些。
五、结论
本文介绍了常微分方程初值问题的两种解法:单步法(Rung-Kula方法)和多步法(Adams方法)。
在实际问题中,如果要求解的精度高,常用的是四阶
Rung-Kula方法,这是因为四阶Rung-Kula方法精度高,编程简单,易于调节步
f x y的次数比同阶的线长,计算过程稳定;其缺点是计算量较大,计算函数(,)
f x y的光滑性较差,Rung-Kula 性多步法(四阶Adams方法)多几倍;另外,如果(,)
f x y较复杂,用Adams 方法的精度还不如改进的Euler方法高。
一般地,若(,)
方法比较合适。
不论采用哪种方法,选取的步长应使落在绝对稳定区域内。
一般在保证精度的条件下,步长尽可能大些,这样可以节省计算量。