常微分方程初值问题数值解法的比较
常微分方程基本概念
常微分方程基本概念常微分方程（Ordinary Differential Equations，简称ODE）是数学分析中的一个重要分支，研究的是一元函数的导数与自变量之间的关系。
它在物理学、工程学、生物学等领域具有广泛的应用。
本文将介绍常微分方程的基本概念和相关知识。
一、常微分方程的定义常微分方程是描述未知函数的导数与自变量之间关系的方程。
一般形式可以表示为：dy/dx = f(x, y)其中，y是未知函数，x是自变量，f(x, y)是已知函数。
二、常微分方程的阶数常微分方程根据未知函数的最高阶导数的阶数不同，可以分为一阶、二阶、高阶等不同阶数的微分方程。
1. 一阶微分方程一阶微分方程是指含有一阶导数的方程。
一般形式可以表示为：dy/dx = f(x, y)例如，y' = 2x + 1就是一个一阶微分方程，其中y'表示y对x的一阶导数。
2. 二阶微分方程二阶微分方程是指含有二阶导数的方程。
一般形式可以表示为：d²y/dx² = f(x, y, dy/dx)例如，y'' + y = 0就是一个二阶微分方程，其中y''表示y对x的二阶导数。
三、常微分方程的初值问题和边值问题常微分方程除了描述函数的导数与自变量之间的关系外，还可以给出一些初始条件或边界条件，从而确定唯一的解。
1. 初值问题初值问题是指在微分方程中给出了函数在某一点的初值条件，要求求解出满足该条件的解。
一般形式可以表示为：dy/dx = f(x, y)，y(x₀) = y₀其中，y(x₀) = y₀表示在点(x₀, y₀)处给定了函数的初始值条件。
2. 边值问题边值问题是指在微分方程中给出了函数在多个点的边界条件，要求求解出满足这些条件的解。
一般形式可以表示为：dy/dx = f(x, y)，y(a) = y_a，y(b) = y_b其中，y(a) = y_a和y(b) = y_b表示在点(a, y_a)和(b, y_b)处给定了函数的边界条件。
浅谈常微分方程初值问题数值解法
浅谈常微分方程初值问题数值解法在自然科学、工程技术、甚至社会科学的一些领域中，常常会遇见一阶常微分方程的求解问题：()上述问题，寻求解的具体表达式十分困难，仅对一些特殊形式的才有可能找到解的解析表达式，在大多情况下，初值问题的解不能用初等函数表示出来即使可写出解的解析表达式，但因为这些表达式过于复杂，要计算它在某些点上的函数值也异常困难。
在实际问题中，经常需要的恰是解在某些点上的函数值，因此研究初值问题的数值解法十分必要。
1 常微分方程初值问题的数值解法常微分方程的近似解法大体可分成三大类：一类是图解法和器械法；第二类是解的近似法；第三类是数值解法，即通过离散化的方法直接求出函数在某些点上的近似值，此数值解仅为精确解的近似解。
其基本原理为：一阶常微分方程的初值问题的解是上变量的连续函数，因此求上述问题的数值解，就是在区间上的若干离散点上用离散化的方法将初值问题化成离散变量的相应问题，从而相应问题的解可作为初值问题理论解的近似值。
由常微分方程的理论可知，只要在区域内连续，且关于满足林普希兹条件，则方程的解存在且唯一。
初值问题的数值解法通常采取“步进法”，而“步进法”又可分为“单步法”和“多步法”两类。
（1）单步法。
所谓“单步法”是指在计算时，只用到前一步的有关信息。
其一般形式为：，主要包括下面三种方法：Euler方法，改进的Euler公式-梯形公式和Runge-Kutta法。
（2）线性多步法。
单步法没有用到前几步计算得到的信息，因此为了提高精度，需重新计算多个点处的函数数值，如RK方法，故计算量较大。
线性多步法的基本思想是充分利用前面的已知信息来构造精度高且计算量小的算法来计算。
多步法常用方法是线性多步法，求解公式为：构造的常用方法是Taylor展开和数值积分方法。
常用的线性多步公式有：四阶Adams显式公式：四阶Adams隐式公式：四阶Milne显式公式：三阶Hamming公式：（隐式公式）预测校正系统和预测校正修正法：一般地，同阶的隐式法比显式法精确，而且数值稳定性好，但隐式公式中的求解较难，需要用到迭代法，这就增加了计算量。
数值分析 第9章 常微分方程初值问题数值解法
9 .2 .2 梯形方法/* trapezoid formula */— 显、隐式两种算法的平均 为得到比欧拉法精度高的计算公式, 在等式( 2.4) 右端积分 中若用梯形求积公式近似, 并用yn 代替y ( xn ) , yn+1 代替y ( xn+1 ) ,则得
h yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn 1 )], 2
yn 1 yn f ( xn , yn ), xn 1 xn
即 yn+1 = yn + hf ( xn , yn ) . ( 2 .1)
这就是著名的欧拉( Euler ) 公式.
• 若初值y0 已知, 则依公式( 2.1)可逐步算出
• y1 = y0 + hf ( x0 , y0 ) ,
为了分析迭代过程的收敛性, 将( 2. 7) 式与(2. 8 )式相减, 得
h ( k 1) (k ) yn 1 yn [ f ( x , y ) f ( x , y 1 n 1 n 1 n 1 n 1 )] 2
于是有
| yn 1 y
( k 1) n 1
hL (k ) | | yn 1 yn 1 |, 2
| f ( x, y1 ) f ( x, y2 ) | L | y1 y2 |, y1, y2 R,
定理1 设f在区域Ｄ＝{(x,y)|a≤x ≤b,y∈R}上连续， 关于y满足利普希茨条件，则对任意x0 ∈[a,b], y0 ∈R，常微分方程初值问题(1.1)式和(1.２)式当x ∈[a,b]时存在唯一的连续可微解y(x). 定理２ 设f在区域Ｄ上连续，且关于y满足利 普希茨条件，设初值问题
1 2 1 2 dy x ydy xdx y x c 2 2 dx y y (0) 2 y2 x2 4
第9章 常微分方程初值问题数值解法
数值分析
第9章 常微分方程初值问题数值解法
《常微分方程》中介绍的微分方程主要有:
(1)变量可分离的方程 (2)一阶线性微分方程(贝努利方程) (3)可降阶的一类高阶方程 (4)二阶常系数齐次微分方程 (5)二阶常系数非齐次微分方程 (6)全微分方程 本章主要介绍一阶常微分方程初值问题的数值解法。
进一步: 令
y n1 y n
xn 1 xn
y n 1 y( x n 1 ) , y n y( x n )
f ( x , y( x ))dx h f ( x n , y n )
宽
9
高
实际上是矩形法
数值分析
第9章 常微分方程初值问题数值解法
(3)
用Taylor多项式近似并可估计误差
解决方法：有的可化为显格式，但有的不行 18
数值分析
第9章 常微分方程初值问题数值解法
与Euler法结合，形成迭代算法 ，对n 0，2， 1，
( yn0 )1 yn hf x n , yn ( k 1) h ( yn1 yn f x n , yn f x n1 , ynk )1 2
7
数值分析
第9章 常微分方程初值问题数值解法
建立数值解法的常用方法
建立微分方程数值解法,首先要将微分方程离散 化. 一般采用以下几种方法: (1) 用差商近似导数
dy yx yx x x dx x y
n 1 n n 1 n
n
,
n
进一步: 令
yn1 y( xn1 ) , yn y( xn )
由 x0 , y0 出发取解曲线 y y x 的切线（存在！），则斜率
常微分方程初值问题数值解法
0.4 1.3582 1.3416 0.9 1.7178 1.6733
0.5 1.4351 1.4142 1.0 1.7848 1.7321
7
初值问题（2.2）有解 y ，1按2这x 个解析式子
算出的准确值 y(x同n )近似值 一y起n 列在表9-1中，两者 相比较可以看出欧拉方法的精度很差.
17
所以，局部截断误差可理解为用方法（2.10）计算一步的 误差，也即公式（2.10）中用准确解y(x代) 替数值解产生
的公式误差.
根据定义，显然欧拉法的局部截断误差
Tn1 y( xn1) y( xn ) hf ( xn , y( xn ))
y(xn h) y(xn ) hy(xn )
y(2) n1
yn
hf
( xn1,
y (1) n1
).
11
如此反复进行，得
y (k 1) n1
yn
hf
( xn1,
y(k) n1
),
(k 0,1, ).
（2.6）
由于 f (x,对y) 满足y 利普希茨条件（1.3）. 由（2.6）减 （2.5）得
y (k 1) n 1
yn1
h
f
( xn1,
y(k) n 1
积分曲线上一点 (x的, y切)线斜率等于函数 值.
的f (x, y)
如果按函数 f (在x, y) 平x面y上建立一个方向场，那 么，积分曲线上每一点的切线方向均与方向场在该点的方 向相一致.
基于上述几何解释，从初始点 P0 (x出0 ,发y0，) 先依 方向场在该点的方向推进到 x 上x1一点 ，P然1 后再从 P1 依方向场的方向推进到 x 上x2一点 ，循P2此前进做出
数学物理方程的数值解法
数学物理方程的数值解法数学物理方程是自然界和科学中描述物体运动、能量转化和相互作用的基本规律。
我们通常使用数值解法来求解这些方程，以得到近似的解析解。
数值解法既可以用于数学问题，也可以用于物理问题。
本文将介绍几种常见的数学物理方程的数值解法。
一、微分方程的数值解法微分方程是描述物体运动和变化的重要工具。
常见的微分方程有常微分方程和偏微分方程。
常见的数值解法包括：1. 欧拉法（Euler's method）欧拉法是最简单的数值解法之一，通过将微分方程离散化为差分方程，在每个小时间步长上近似计算微分方程的导数。
欧拉法易于实现，但精度相对较低。
2. 龙格-库塔法（Runge-Kutta method）龙格-库塔法是一类常用的数值解法，包括二阶、四阶等不同的步长控制方法。
龙格-库塔法通过计算多个离散点上的导数来近似微分方程，精度较高。
3. 有限差分法（Finite difference method）有限差分法是一种常用的数值解法，将微分方程转化为差分方程并在网格上逼近微分方程的导数。
有限差分法适用于边值问题和初值问题，且精度较高。
二、积分方程的数值解法积分方程描述了给定函数的积分和积分变换之间的关系。
常见的数值解法有：1. 数值积分法数值积分法是通过数值逼近求解积分方程，常用的数值积分法包括梯形法则、辛普森法则等。
数值积分法适用于求解一维和多维积分方程。
2. 蒙特卡洛法（Monte Carlo method）蒙特卡洛法通过随机采样和统计分析的方法，将积分方程转化为概率问题，并通过大量的随机样本来估计积分值。
蒙特卡洛法适用于高维空间和复杂积分方程。
三、优化问题的数值解法优化问题是寻找在给定约束条件下使目标函数取得极值的数学问题。
常见的数值解法有：1. 梯度下降法（Gradient descent method）梯度下降法是一种常用的优化算法，通过迭代和梯度方向来寻找目标函数的局部最优解。
梯度下降法适用于连续可导的优化问题。
微分方程的一些通解和初值问题的解法
微分方程的一些通解和初值问题的解法微分方程作为数学中一个极其重要的分支，它具有广泛的应用背景，包括自然科学、工程技术等多个领域中都有着广泛的应用。
微分方程的求解则是这门学科中一个很关键的问题，尤其是对于一些实际问题，其初值条件决定了微分方程的具体解，本文将探讨一些微分方程的通解以及初值问题解法。
1. 常微分方程的通解对于一个n阶常微分方程，如果它可以表示为：$$F\Bigg(x,\frac{dy}{dx},\frac{d^2 y}{dx^2},\cdots,\frac{d^ny}{dx^n}\Bigg)=0$$其中$y$是自变量$x$的函数，则这个方程是一个n阶常微分方程。
对于这类方程，可以根据它的阶数以及特点进行分类求解。
(1)一阶常微分方程通解这类方程形式如下：$$\frac{dy}{dx}=f(x,y)$$其中$f(x,y)$是定义在某个区域上的函数。
对于这类方程，我们可以通过分离变量的方式进行求解，即：$$\frac{dy}{f(x,y)}=dx$$两边同时积分得到：$$\int\frac{1}{f(x,y)}dy=\int dx+C$$其中$C$是积分常数，通过这个式子可以求得$y$的通解。
(2)二阶常微分方程通解这类方程形式如下：$$y''+p(x)y'+q(x)y=f(x)$$其特点是含有二阶导数项，可用特征方程进行求解。
将一般形式二阶常微分方程的通解表示为$y=c_1y_1+c_2y_2$，其中$c_1$和$c_2$是常数，$y_1$和$y_2$是方程的解，满足$y_1$和$y_2$的任意线性组合都是方程的解。
如果解$y_1$和$y_2$线性无关，则它们构成了二阶常微分方程的通解。
(3)n阶常微分方程通解通常情况下，n阶常微分方程表示为：$$y^{(n)}+a_{n-1}(x)y^{(n-1)}+\cdots+a_1(x)y'+a_0(x)y=f(x)$$我们可以通过求解$n$次的导数，得到这个方程的通解。
常微分方程初值问题的数值解法
常微分方程初值问题数值解法初值问题：即满足初值条件的常微分方程的解y′=f(x,y),x∈[x0,b]y(x0)=y0.定理1（利普希茨条件）若存在正数L，使得对任意，y1，y2，有|f(x,y1)−f(x,y2)|≤L|(y1−y2)|定理2（解存在性）①若函数f在方区域x∈[a,b],y∈R连续，②函数f关于y 满足利普希茨条件，则对任意x∈[a,b]，常微分方程存在唯一的连续可微数值解.两类问题：①单步法---计算下一个点的值yn+1只需要用到前面一个点的值yn②多步法---计算下一个点的值yn+1需要用到前面l个点的值yl1、欧拉法---下一个点的计算值等于前一个点的计算值加上步长乘以前一个点的函数值•具体过程一些批注：显式欧拉方程指下一步要计算的值，不在迭代方程中；隐式欧拉方程指下一步要计算的值，在迭代方程中。
怎么计算隐式欧拉方程----要借助显示欧拉迭代计算---一般用迭代法-----迭代---将微分方程在区间[xn,xn+1]进行积分，然后函数f进行近似，即可得到迭代方程-----迭代方程收敛性？由函数关于y满足利普希茨条件，可以推出迭代公式收敛。
•局部截断误差：假设前n步误差为0，我们计算第n+1步的误差，将次误差称为局部截断误差，且局部误差为O(hp+1)•p阶精度：由理论证明：若局部误差阶的时间复杂度为O(hp+1)，则整体误差阶为O(hp)我们称公式精度为p。
•显示欧拉法与隐式欧拉法•梯形方法----将显式欧拉迭代方程与隐式欧拉迭代方程做一下加权平均，构造的计算公式.•改进的欧拉方法---思想：因为梯形公式是隐式公式，将显式欧拉公式对下一步的计算值进行预估，用梯形公式对下一步的计算值进行校正.2、龙格-库塔方法思想：根据Lagrange中值定理，下一次的计算值可以用前一次的计算值加上h乘以前一个点的斜率；而这个斜率用该区间上的多个点的斜率的算数平均来逼近。
注意：怎么计算任意斜率Ki？第i个点的斜率Ki有微分方程可以算出f′=f(xn,yn)所以要算的f(xn,yn)值，由欧拉法即可算出, yn+1=yn+hf′•2阶-龙格-库塔方法----类似改进的欧拉法根据Lagrange中值定理，下一次的计算值可以用前一次的计算值加上h乘以斜率；而这个斜率用区间上的端点和中点的斜率的算数平均来逼近。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值计算实践—课程设计报告
课题名称
常微分方程初值问题数值解法的比较
完成时间
2013-1-17
姓名
班级
学号
成绩
一. 实验目的及内容
1实验目的：（1）了解常微分方程初值问题的理论背景以及初值问题稳定性、收敛性的研究；
（2）熟练掌握欧拉法、改进欧拉法、龙格-库塔法以及截断误差分析；
输入：[x,y]=Gaijineuler('f',0,1,1,0.1)
结果：x =
0
0.1000
0.2000
0.3000
0.4000
0.5000
0.6000
0.7000
0.8000
0.9000
1.0000
y =
1.0000
1.0959
1.1841
1.2662
1.3434
1.4164
1.4860
1.5525
end
R=[T' Y'];
（四阶龙格-库塔）
functionR=Longgekuta4(f,a,b,aZ,h)
%a，b为端点，h为步长，aZ为初值
n=(b-a)/h;
T=zeros(1,n+1);%定义向量
Y=zeros(1,n+1);
T=a:h:b;%计算各个分点
Y(1)=aZ;%初值赋予
forj=1:n
0.3000 1.3411
0.4000 1.4711
0.5000 1.6082
0.6000 1.7520
0.7000 1.9021
0.8000 2.0580
0.9000 2.2195
1.0000 2.3863
五．计算结果分析:
方法
显示欧拉
简单
精度低
隐式欧拉
稳定性最好
精度低，计算量大
梯形公式
精度提高
计算量大
（4）对右端 用右矩形公式得 ，也叫隐式欧拉法。
误差分析：1.称 为计算 时的局部截断误差；
2.如果数值方法的局部截断误差为 ，那么称这种数值方法的阶数是p，其实p为非负整数。通常情况下，步长h越小，p越高，则局部截断误差越小，计算精
初泰勒展开有
则有 可见欧拉方法是一阶方法，精度不是很高。
2．改进欧拉方法：
考虑一阶常微分方程的初值问题 ，如果存在实数 使得 则称f关于y满足利普希茨条件，L称为利普希茨常数。
对于常微分方程初值问题 ，考虑初值 的扰动是问题的解 发生偏差的情形。若 时 的偏差被控制在有界范围内，则称该初值问题是稳定的，否则该初值问题不稳定的。
特别地，若 时 的偏差收敛于零，则称该初值问题是渐进稳定的。
对于初值问题 稳定性的研究，易知其准确解为 ，假定初值经过扰动后变为 ，对于扰动后的解为 因此带来的扰动误差为 ，因此考虑 时 的值，它取决于 。易知，若 ，则原问题是稳定的；若 ，则原问题是不稳定的；若 ，则原问题是渐进稳定的。
实际遇到的大多数常微分方程初值问题都是稳定的，因此在后面的讨论数值解法时这常常是默认条件。
y(i+1)=(yp+yc)/2;%根据改进欧拉公式计算结果
end
x=x';
y=y';
1．
3.龙格-库塔代码如下：
（三阶龙格-库塔）
functionR=Longgekuta3(f,a,b,aZ,h)
%a，b为端点，h为步长，aZ为初值
n=(b-a)/h;
T=zeros(1,n+1);%定义向量
Y=zeros(1,n+1);
2.稳定性：
例：考察初值问题 在区间[0,0.5]上的解，分别用欧拉显，隐式格式，改进欧拉格式计算数值解：
节点
欧拉显式
欧拉隐式
改进欧拉法
精确解
0
1.0000
1.0000
1.0000
1.0000
0.1
-2.0000
2.5*10^(-1)
2.5000
4.9787*10^(-2)
0.2
4.0000
6.2500*10^(-2)
T=a:h:b;%计算各个分点
Y(1)=aZ;%初值赋予
forj=1:n
k1=feval(f,T(j),Y(j));
k2=feval(f,T(j)+h/2,Y(j)+k1*h/2);
k3=feval(f,T(j)+h,Y(j)-h*k1+k2*2*h);
Y(j+1)=Y(j)+(k1+4*k2+k3)*h/6;%根据公式计算
中点公式
精度提高，显式
多一个初值，可能影响精度
1.收敛性：
若某算法对于任意固定的x=xi=x0+i h，当h0 (同时i)时有yiy(xi)，则称该算法是收敛的。
以下讨论的都是单步法（指在计算 时只用到它前一步的信息 ）欧拉法，龙格-库塔法都是单步法的例子；
例：就初值问题 考察欧拉显式格式的收敛性。
二．相关背景知识介绍以及初值问题稳定性的研究：
在科学与工程问题中，常微分方程表述物理量的变化规律，应用非常广泛，比如，天体运动的轨迹，机器人控制，化学反应过程的描述和控制以及电路瞬态过程分析等。这些问题中要求解随时间变化的物理量，即位置函数 表示时间，而微分方程描述了未知函数与它的一阶或高阶导数之间的关系。
6.2500
2.4788*10^(-3)
0.3
-8.0000
1.5625*10^(-2)
1.5625*10^(1)
1.2341*10^(-4)
0.4
1.6000*10^(1)
3.9063*10^(-3)
3.9063*10^(1)
6.1442*10^(-6)
0.5
-3.2000*10^(1)
9.7656*10^(-1)
0.7000 1.9487
0.8000 2.1436
0.9000 2.3579
1.0000 2.5937
ans =
Columns 1 through 7
1.0000 1.1000 1.2100 1.3310 1.4641 1.6105 1.7716
2.1436 2.3579 2.5937
1.6165
1.6782
1.7379
龙格-库塔计算结果：
（三阶）输入：Longgekuta3('ff',0,1,1,0.1)
结果：ans =
Longgekuta3('ff',0,1,1,0.1)
ans =
0 1.0000
0.1000 1.1048
0.2000 1.2188
0.3000 1.3411
0.4000 1.4711
end
R=[T' Y'];
四．数值结果：
输入：定义M文件并保存ffx.m
Euler1('ffx',0,1,1,10)
结果：T =
0 1.0000
0.1000 1.1000
0.2000 1.2100
0.3000 1.3310
0.4000 1.4641
0.5000 1.6105
0.6000 1.7716
梯形公式：对右端 用梯形公式得+显然梯形公式是隐式公式。
改进欧拉公式：先用欧拉公式求的一个初步的近似值 ，成为预测值，预测值 的精度可能达不到要求，在用梯形公式将他校正一次，记为 ，这个结果成为校正值。
预测：
校正：
误差分析：记 为改进欧拉公式在 处的截断误差，
记 因此有
， 表示在 出的局部截断误差。由此得，梯形公式的局部截断误差为 ，因此改进欧拉的截断误差为 ，可见改进欧拉的方法是二阶方法，改进欧拉方法优于欧拉方法。
k1=feval(f,T(j),Y(j));
k2=feval(f,T(j)+h/2,Y(j)+k1*h/2);
k3=feval(f,T(j)+h/2,Y(j)+k2*h/2);
k4=feval(f,T(j)+h,Y(j)+k3*h);
Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)*h/6;%根据公式计算
向前欧拉公式结果：
输入：Euler2('ffx',0,1,1,10)
结果T =
0 1.0000
0.1000 1.1110
0.2000 1.2344
0.3000 1.3716
0.4000 1.5240
0.5000 1.6933
0.6000 1.8814
0.7000 2.0904
0.8000 2.3227
3．龙格—库塔法:
根据拉格朗日微分中值定理， ，记 得到 ，这样，只要给出一种计算 的算法，就可以得到相应的计算公式。欧拉公式可以写为
改进欧拉公式可以写成 因此推出一般的推出广式
称为p阶龙格-库塔方法，简称p阶R-K方法。
因为 这里的 均为常数。
因为给定的系数不唯一，因此这里的常数有无穷多个解，下面是特殊情况下和一般情况下得结果
（一阶龙格-库塔）当r=1时，这就是欧拉法。
（二阶龙格-库塔）当r=2时， ，这就是改进欧拉法。
三阶和四阶龙格-库塔也只是在一般情况下得结果。
三阶
四阶
其局部截断误差是：
三．程序代码
欧拉法代码如下：
（1）向前欧拉法：
functiony=Euler1(fun,x0,y0,xN,N)
%fun为一阶微分方程，x0，y0为初始条件，xN为取值范围的一个端点，N为区间个数