常微分方程数值解
- 格式：ppt
- 大小：1.34 MB
- 文档页数：24
下载文档原格式
合集下载
数值分析常微分方程数值解法
7
第8页/共105页
➢ 数值积分方法(Euler公式)
设将方程 y=f (x, y)的两端从 xn 到xn+1 求积分, 得
y( xn1) y( xn )
xn1 f ( x, y( x))dx :
xn
xn1 F ( x)dx
xn
用不同的数值积分方法近似上式右端积分, 可以得到计算 y(xn+1)的不同的差分格 式.
h2 2
y''( )
Rn1
:
y( xn1)
yn1
h2 2
y''( )
h2 2
y''( xn ) O(h3 ).
局部截断误差主项
19
第20页/共105页
➢ 向后Euler法的局部截断误差
向后Euler法的计算公式
yn1 yn hf ( xn1, yn1 ), n 0, 1, 2,
定义其局部截断误差为
y 计算 的n递1 推公式，此类计算格式统称为差分格式.
3
第4页/共105页
数值求解一阶常微分方程初值问题
y' f ( x, y), a x b,
y(a)
y0
难点: 如何离散 y ?
➢ 常见离散方法
差商近似导数 数值积分方法 Taylor展开方法
4
第5页/共105页
➢ 差商近似导数(Euler公式)
(0 x 1)
y(0) 1.
解 计算公式为
yn1
yn
hfn
yn
h( yn
2xn ), yn
y0 1.0
n 0, 1, 2,
取步长h=0.1, 计算结果见下表
13
第8页/共105页
➢ 数值积分方法(Euler公式)
设将方程 y=f (x, y)的两端从 xn 到xn+1 求积分, 得
y( xn1) y( xn )
xn1 f ( x, y( x))dx :
xn
xn1 F ( x)dx
xn
用不同的数值积分方法近似上式右端积分, 可以得到计算 y(xn+1)的不同的差分格 式.
h2 2
y''( )
Rn1
:
y( xn1)
yn1
h2 2
y''( )
h2 2
y''( xn ) O(h3 ).
局部截断误差主项
19
第20页/共105页
➢ 向后Euler法的局部截断误差
向后Euler法的计算公式
yn1 yn hf ( xn1, yn1 ), n 0, 1, 2,
定义其局部截断误差为
y 计算 的n递1 推公式，此类计算格式统称为差分格式.
3
第4页/共105页
数值求解一阶常微分方程初值问题
y' f ( x, y), a x b,
y(a)
y0
难点: 如何离散 y ?
➢ 常见离散方法
差商近似导数 数值积分方法 Taylor展开方法
4
第5页/共105页
➢ 差商近似导数(Euler公式)
(0 x 1)
y(0) 1.
解 计算公式为
yn1
yn
hfn
yn
h( yn
2xn ), yn
y0 1.0
n 0, 1, 2,
取步长h=0.1, 计算结果见下表
13
第7章 常微分方程数值解法
来描述。式中，x是质量m离平衡位置（0点）的距离；
2
t 是时间 ；c是弹簧常数。
3
当弹簧在振动开始时刻t = t0 时的初始位置x( t0)= x0 和初速度
dx dt
x (t ) x
0 0 t t0
确定时，弹簧的振动规律x(t) 也就唯一确定。这就 是一个常微分方程的初值问题，可写成：
称为步长。当
h h （常值）时称为等步长，有
i
xi x0 ih, (i 1,2,...n)
或
xi1 xi h, (i 0,1,2,...n 1)
7
因为初值问题中的初始条件 y (a) y0 为已知，即 可利用已知的 y 0 来求出下一节点处 y ( x1 ) 的近似值 y 1
初值问题（1.1）的数值解法，常采用差分方法，即把
一个连续的初值问题离散化为一个差分方程来求解。即将 （1.1）离散化后，求找其解y
= y (x)在一系列离散节点
6
a = x0 < x1 < … < xi < … < xn = b
上的近似值y0, y1, …, yn。
两相邻节点间的距离
hi = xi+1 - xi (i=0,1,2,…,n-1)
yi 1 yi hf ( xi , y i )(i 0,1,2 , n 1) (2.2)
这就是欧拉（Euler）公式，又称欧拉格式。利用它 可由已知的初值 y0 出发，逐步算出 y1 , y 2 , y n 。这 类形式的方程也称为差分方程。 当假定 yi 为准确，即在 yi y ( xi ) 的前提下来估 计误差 y ( xi 1) y i 1 ，这种截断误差称为局部截断误 差。由（2.1）和（2.2）可知，欧拉格式在节点x i 1 处的局部截断误差显然为：
2
t 是时间 ；c是弹簧常数。
3
当弹簧在振动开始时刻t = t0 时的初始位置x( t0)= x0 和初速度
dx dt
x (t ) x
0 0 t t0
确定时，弹簧的振动规律x(t) 也就唯一确定。这就 是一个常微分方程的初值问题，可写成：
称为步长。当
h h （常值）时称为等步长，有
i
xi x0 ih, (i 1,2,...n)
或
xi1 xi h, (i 0,1,2,...n 1)
7
因为初值问题中的初始条件 y (a) y0 为已知，即 可利用已知的 y 0 来求出下一节点处 y ( x1 ) 的近似值 y 1
初值问题（1.1）的数值解法，常采用差分方法，即把
一个连续的初值问题离散化为一个差分方程来求解。即将 （1.1）离散化后，求找其解y
= y (x)在一系列离散节点
6
a = x0 < x1 < … < xi < … < xn = b
上的近似值y0, y1, …, yn。
两相邻节点间的距离
hi = xi+1 - xi (i=0,1,2,…,n-1)
yi 1 yi hf ( xi , y i )(i 0,1,2 , n 1) (2.2)
这就是欧拉（Euler）公式，又称欧拉格式。利用它 可由已知的初值 y0 出发，逐步算出 y1 , y 2 , y n 。这 类形式的方程也称为差分方程。 当假定 yi 为准确，即在 yi y ( xi ) 的前提下来估 计误差 y ( xi 1) y i 1 ，这种截断误差称为局部截断误 差。由（2.1）和（2.2）可知，欧拉格式在节点x i 1 处的局部截断误差显然为：
常微分方程初值问题数值解法
根据微分方程的性质和初始条件，常 微分方程初值问题可以分为多种类型， 如一阶、高阶、线性、非线性等。
数值解法的必要性
实际应用需求
许多实际问题需要求解常微分方程初值问题，如物理、 化学、生物、工程等领域。
解析解的局限性
对于复杂问题，解析解难以求得或不存在，因此需要 采用数值方法近似求解。
数值解法的优势
未来发展的方向与挑战
高精度算法
研究和发展更高精度的算法，以提高数值解的准确性和稳定性。
并行计算
利用并行计算技术，提高计算效率，处理大规模问题。
自适应方法
研究自适应算法，根据问题特性自动调整计算精度和步长。
计算机技术的发展对数值解法的影响
1 2
硬件升级
计算机硬件的升级为数值解法提供了更强大的计 算能力。
它首先使用预估方法（如欧拉方法）得到一个 初步解，然后使用校正方法（如龙格-库塔方法） 对初步解进行修正，以提高精度。
预估校正方法的优点是精度较高，且计算量相 对较小，适用于各种复杂问题。
步长与误差控制
01
在离散化过程中，步长是一个重要的参数，它决定 了离散化的精度和计算量。
02
误差控制是数值逼近的一个重要环节，它通过设定 误差阈值来控制计算的精度和稳定性。
能够给出近似解的近似值，方便快捷，适用范围广。
数值解法的历史与发展
早期发展
早在17世纪，科学家就开始尝 试用数值方法求解常微分方程。
重要进展
随着计算机技术的发展，数值 解法在20世纪取得了重要进展， 如欧拉法、龙格-库塔法等。
当前研究热点
目前，常微分方程初值问题的 数值解法仍有许多研究热点和 挑战，如高精度算法、并行计
软件优化
软件技术的发展为数值解法提供了更多的优化手 段和工具。
数值解法的必要性
实际应用需求
许多实际问题需要求解常微分方程初值问题，如物理、 化学、生物、工程等领域。
解析解的局限性
对于复杂问题，解析解难以求得或不存在，因此需要 采用数值方法近似求解。
数值解法的优势
未来发展的方向与挑战
高精度算法
研究和发展更高精度的算法，以提高数值解的准确性和稳定性。
并行计算
利用并行计算技术，提高计算效率，处理大规模问题。
自适应方法
研究自适应算法，根据问题特性自动调整计算精度和步长。
计算机技术的发展对数值解法的影响
1 2
硬件升级
计算机硬件的升级为数值解法提供了更强大的计 算能力。
它首先使用预估方法（如欧拉方法）得到一个 初步解，然后使用校正方法（如龙格-库塔方法） 对初步解进行修正，以提高精度。
预估校正方法的优点是精度较高，且计算量相 对较小，适用于各种复杂问题。
步长与误差控制
01
在离散化过程中，步长是一个重要的参数，它决定 了离散化的精度和计算量。
02
误差控制是数值逼近的一个重要环节，它通过设定 误差阈值来控制计算的精度和稳定性。
能够给出近似解的近似值，方便快捷，适用范围广。
数值解法的历史与发展
早期发展
早在17世纪，科学家就开始尝 试用数值方法求解常微分方程。
重要进展
随着计算机技术的发展，数值 解法在20世纪取得了重要进展， 如欧拉法、龙格-库塔法等。
当前研究热点
目前，常微分方程初值问题的 数值解法仍有许多研究热点和 挑战，如高精度算法、并行计
软件优化
软件技术的发展为数值解法提供了更多的优化手 段和工具。
常微分方程的解析解与数值解
常微分方程的解析解与数值解常微分方程是数学中的重要分支，广泛应用于物理学、工程学、生物学等领域。
解析解和数值解是求解常微分方程的两种常用方法。
本文将介绍常微分方程的解析解和数值解的概念、特点以及应用，并讨论两者在不同情况下的优缺点。
一、解析解解析解是指通过数学方法直接获得的方程的解。
对于某些简单的常微分方程，可以通过变量分离、分离变量、常数变易等方法获得解析解。
解析解具有以下几个特点：1. 精确性：解析解是通过数学方法得到的，是方程的精确解。
它可以给出方程在任意时刻的解，对于问题的研究具有重要意义。
2. 通用性：解析解适用于一类具有相同形式的常微分方程。
一旦求得了一类方程的解析解，就可以应用到同类问题中。
3. 物理含义明确：解析解通常具有明确的物理含义，能够帮助我们理解问题的本质和规律。
解析解在一些特定情况下具有明显的优势。
例如，当方程形式简单、边界条件明确、初值明确时，解析解能够提供准确的结果。
此外，解析解也有助于我们对问题的理论分析和深入研究。
二、数值解数值解是通过数值计算方法获得的方程的近似解。
对于复杂的常微分方程，往往难以找到解析解，这时候就需要借助数值方法进行求解。
数值解具有以下几个特点：1. 近似性：数值解是通过数值计算获得的，只能提供方程的近似解。
随着计算步长的减小，近似解会逐渐接近真实解。
2. 条件灵活：数值解对问题的条件要求相对较低。
例如，数值方法可以求解非线性方程、高阶方程、边值问题等各种复杂情况。
3. 计算复杂度：数值解通常需要借助计算机进行迭代计算，计算复杂度较高。
数值解在实际问题中应用广泛且有效。
数值方法可以通过逼近、插值、差分等数值计算技术，将方程转化成逐步计算的步骤，获得精确度可控的近似解。
数值解的优势在于对于复杂问题的求解能力和计算相对高效。
三、解析解与数值解的比较解析解和数值解各自具有不同的特点和优势，在不同的问题和求解需求中有着相应的应用。
解析解在以下情况下具有优势：1. 简单线性方程：对于形式简单的一阶线性常微分方程，如首次线性方程、可分离变量方程等，可以通过解析方法求得解析解。
常微分方程数值解法2线性多步法
对于线性多步法，其收敛性取决于微分方程的解的性质和方法的阶数。一般来说，高阶方法具有更好 的收敛性。
03
常见的线性多步法
欧拉方法
总结词
欧拉方法是常微分方程数值解法中最简单的一种方法，它基于线性近似，通过已知的函 数值来估计新的函数值。
详细描述
欧拉方法的基本思想是利用已知的函数值来估计下一个点的函数值。具体来说，假设我 们有一个函数 (y = f(x))，在已知 (x_0) 处的函数值 (y_0 = f(x_0)) 的情况下，欧拉方法 通过线性插值来估计 (x_1) 处的函数值 (y_1)，即 (y_1 = y_0 + h f(x_0))，其中 (h) 是
05
线性多步法的优缺点
优点
稳定性好
线性多步法在处理常微分方程时具有较好的数值稳定性， 能够有效地抑制数值振荡，提高计算结果的精度。
01
易于实现
线性多步法的计算过程相对简单，易于 编程实现，适合于大规模数值计算。
02
03
精度可调
通过选择不同的步长和线性多步法公 式，可以灵活地调整计算结果的精度， 满足不同的数值模拟需求。
改进方法的收敛性
研究收敛性条件
深入研究线性多步法的收敛性条件，了解哪些情况下方法可能不收 敛，并寻找改进措施。
优化迭代算法
通过优化迭代算法，提高方法的收敛速度和精度，减少迭代次数， 提高计算效率。
引入预处理技术
利用预处理技术对线性系统进行预处理，改善系统的条件数，提高方 法的收敛性。
拓展应用领域
在工程问题中的应用
控制系统设计
在工程领域中，线性多步法可以用于控制系统设计，通过 建立控制系统的数学模型，设计控制算法和控制器，实现 系统的稳定性和性能优化。
03
常见的线性多步法
欧拉方法
总结词
欧拉方法是常微分方程数值解法中最简单的一种方法，它基于线性近似，通过已知的函 数值来估计新的函数值。
详细描述
欧拉方法的基本思想是利用已知的函数值来估计下一个点的函数值。具体来说，假设我 们有一个函数 (y = f(x))，在已知 (x_0) 处的函数值 (y_0 = f(x_0)) 的情况下，欧拉方法 通过线性插值来估计 (x_1) 处的函数值 (y_1)，即 (y_1 = y_0 + h f(x_0))，其中 (h) 是
05
线性多步法的优缺点
优点
稳定性好
线性多步法在处理常微分方程时具有较好的数值稳定性， 能够有效地抑制数值振荡，提高计算结果的精度。
01
易于实现
线性多步法的计算过程相对简单，易于 编程实现，适合于大规模数值计算。
02
03
精度可调
通过选择不同的步长和线性多步法公 式，可以灵活地调整计算结果的精度， 满足不同的数值模拟需求。
改进方法的收敛性
研究收敛性条件
深入研究线性多步法的收敛性条件，了解哪些情况下方法可能不收 敛，并寻找改进措施。
优化迭代算法
通过优化迭代算法，提高方法的收敛速度和精度，减少迭代次数， 提高计算效率。
引入预处理技术
利用预处理技术对线性系统进行预处理，改善系统的条件数，提高方 法的收敛性。
拓展应用领域
在工程问题中的应用
控制系统设计
在工程领域中，线性多步法可以用于控制系统设计，通过 建立控制系统的数学模型，设计控制算法和控制器，实现 系统的稳定性和性能优化。
常微分方程数值解法
第七章 常微分方程数值解法常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。
另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。
因此，有必要探讨常微分方程初值问题的数值解法。
本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。
第一节 欧拉法求解常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy（1）的数值解，就是寻求准确解)(x y 在一系列离散节点 <<<<<n x x x x 210 上的近似值 ,,,,,210n y y y y{}n y 称为问题的数值解，数值解所满足的离散方程统称为差分格式，1--=i i ix x h 称为步长，实用中常取定步长。
显然，只有当初值问题(1)的解存在且唯一时，使用数值解法才有意义，这一前提条件由下 面定理保证。
定理 设函数()y x f ,在区域+∞≤≤-∞≤≤y b x a D ,:上连续，且在区域D 内满足李普希兹(Lipschitz)条件，即存在正数L ，使得对于R 内任意两点()1,y x 与()2,y x ,恒有()()2121,,y y L y x f y x f -≤-则初值问题(1)的解()x y 存在并且唯一。
一、欧拉法(欧拉折线法)若将函数)xy 在点nx处的导数()n x y '用两点式代替， 即()hx y x y x y n n n )()(1-≈'+，再用n y 近似地代替()n x y ,则初值问题(1)变为⎩⎨⎧==++=+ ,2,1,0),()(001n x y y y x hf y y n n n n(2)(2)式就是著名的欧拉(Euler)公式。
以上方法称 为欧 拉法或欧拉折线法。
《常微分方程的数值解法》论文
《常微分方程的数值解法》论文《常微分方程的数值解法》常微分方程（ODE）是研究物理过程的重要工具，其伴随着极大的应用价值。
当一个物理系统被简化为一个常微分方程，它就可以用于描述物理学中的各种现象。
但是，大多数现实系统的常微分方程未能得到解析解，因此，数值解法就变得非常重要。
本文将研究并比较几种常见的常微分方程数值解法，诸如Euler法、奇异点法、Runge-Kutta法、前向差分法等，以便更好地提供协助解决常微分方程。
首先，Euler法是常用的数值解法之一，它主要用于解决常微分方程模型。
其核心思想是将微分方程通过采用不断变化的步长对状态量求近似值，并通过预测下一步的值来求解微分方程，从而达到求解常微分方程的目的，且操作简单、容易理解。
但是，由于其步长的不动性，往往使得其精度较低，因此，当遇到复杂环境时，Euler法的表现就有些不尽如人意。
此外，另一种常见的数值解法是奇异点法。
此法将一个微分方程情况分解成多个分段函数，每一段函数都可以精确求解，从而可以求解复杂的微分方程。
它的特点是分段的每一部分的精度和复杂度都较低，而且运行效率也较快，但是，奇异点法的精度需要在段间合理设定，然后再进行微调，以保证数值模拟的准确性。
其次，Runge-Kutta法是一种常用的数值解法，它可以有效地求解一些常微分方程，其原理是利用积分函数插值，然后利用积分函数求近似值，最后根据边界条件求取解析结果。
Runge-Kutta法的步长可以随着计算过程的进行而逐步变化，这样可以使得误差得到有效控制，而且可以有效地控制误差，保证算法精度，但是由于其计算效率较低，因此在求解复杂的常微分方程时，Runge-Kutta法的表现并不尽人意。
最后，前向差分法是一种求解常微分方程的数值解法，它利用求取未知函数的一阶导数和二阶导数的值，然后通过求解一次和二次中点差分的方式，从而得到数值解。
它的有点是能够得到较高的精确度，且即使步长变化时也可以控制误差，但前向差分法要求在微分方程中必须有高阶导数，这就要求微分方程是复杂的，除此之外，除了必须计算高次导数外，它的计算量也比较大。
第5章常微分方程数值解法
2hyxn
2h3 3!
y
yxn1
yxn1 2hyxn
h3 3
y ②
将①、②两式相减：
y
xn1
h3 yn1 3
y
——两步法局部截断误差
18
第19页/共24页
2024年8月7日
（2）梯形公式
yn1
yn
h 2
f
xn , yn
f
xn1 , yn1
yxn
hf
2
xn ,
yn1 yxn1 2hyxn ①
2024年8月7日
17
第18页/共24页
2024年8月7日
将函数用泰勒级数展开：（ h 较小， 相差不大）
yxn1
yxn
h 1!
yxn
h2 2!
yxn
h3 3!
y
yxn1
yxn
h
1!
yxn
h2
2!
yxn
h3
3!
y
yxn1
yxn1
yxn
f
xn1 , yxn1
yn1
yxn
h 2
y
xn
yxn1 ①
将函数用泰勒级数展开：
yxn1
yxn
h 1!
yxn
h2 2!
yxn
h3 3!
y ②
yxn1
yxn
h 1!
yxn
h2 2!
y ③
（ h 较小， 相差不大）
19
第20页/共24页
2024年8月7日
①、②两式相减，并代入③式：
（图示表示梯形法计算结果）
第5章_常微分方程数值解法
（5.2.6）
由于方程关于 uk +1 是隐式形式，所以式（5.2.6）称为隐式 Euler 公式。前面显式和隐式 Euler 公式在计
u '(tk ) ≈
得到的递推公式：
u (tk +1 ) − u (tk −1 ) 2h
（5.2.7）
uk +1 ≈ uk −1 + 2hf (tk , uk )
在计算 uk +1 时，需要用到前两步结果 uk −1 , uk ，称为两步法公式。 （2）积分近似方法 将（5.2.1）式的微分方程写成 du = f (t , u )dt ，在区间 [tk , tk +1 ] 上积分，有：
5.2.2 Runge-Kutta 方法 Euler 方法比较简单，但它的收敛阶数低。可以利用 Taylor 展开式构造高阶的单步方法。Euler 公式 可以看成是由一阶 Taylor 展开式得到的，所以应用高阶 Taylor 展开就可以得到高阶单步法。例如：将 u (tk +1 ) 在 tk 处作 q 阶 Taylor 展开：
dy = a − by (t ) dt
是一阶常微分方程，而
2 ∂ 2 u ( x, t ) 2 ∂ u ( x, t ) a = ∂t 2 ∂x 2
（5.1.1）
（5.1.2）
是二阶偏微分方程。 所有使微分方程成为等式的函数，都是微分方程的解；在 n 阶微分方程中，将微分方程的含有 n 个任 意常数的解称为该微分方程的通解。为确定微分方程通解中的任意常数而需要的条件称为定解条件；定解 条件可以分为初始条件和边界条件两类。由微分方程和定解条件一起构成的问题称为微分方程定解问题。 根据定解条件的不同，常微分方程分为初值问题和边值问题；若定解条件是描述函数在一点（或初始 点）处状态的，则称为初值问题，一阶常微分方程初值问题的一般形式为：
高等工程数学第五讲,常微分方程数值解法
n h2 M (1 hL)k 2 k 0
h2 M (1 hL)n 1 1 hM (1 hL) n 1 1 2 (1 hL) 1 2 L
17
由 (n 1) h b a (1 hL)n 1 (1 hL)
b a h
而e hL
yn 1
1 yn [ f ( xn , yn ) f ( xn 1 , yn 1 )] 2
梯形公式的局部截断误差 h3 Rn 1 y ( xn 1 ) yn 1 y( ), xn xn 1 12 梯形公式为二阶方法。
10
二, 改进Euler法(梯形法) 由于梯形公式仍然是隐式方法，一般用迭代法 求解，若每一步只代一次，相当将Euler公式与梯 形公式结合使用。 x n ,y n ) 预测 yn 1 yn hf ( (10.9) h yn 1 yn [f (x n , y n ) f (x n , 校正 1y n )] 1 2 为上机编程方便，改写为 y p yn hf ( xn, yn) yq yn hf ( xn h, y p) yn 1 1 ( y p yq ) 2
y x y ,0 x 1 y(0) 0
5
解： 易知精确解为y x e x 1,
Euler公式的具体形式为yn+1 =yn +h(xn -yn ),(n=0,1, )
若取h=0.1,N=10,可由初值y0=0逐次算出.见下表
xn 0 0.1 0.2 0.3 0.9 1.0
f ( x, y) f ( x, y ) L y y ,(10.2)
由常微分方程理论知，初值问题（10.1） 有解且唯一。
h2 M (1 hL)n 1 1 hM (1 hL) n 1 1 2 (1 hL) 1 2 L
17
由 (n 1) h b a (1 hL)n 1 (1 hL)
b a h
而e hL
yn 1
1 yn [ f ( xn , yn ) f ( xn 1 , yn 1 )] 2
梯形公式的局部截断误差 h3 Rn 1 y ( xn 1 ) yn 1 y( ), xn xn 1 12 梯形公式为二阶方法。
10
二, 改进Euler法(梯形法) 由于梯形公式仍然是隐式方法，一般用迭代法 求解，若每一步只代一次，相当将Euler公式与梯 形公式结合使用。 x n ,y n ) 预测 yn 1 yn hf ( (10.9) h yn 1 yn [f (x n , y n ) f (x n , 校正 1y n )] 1 2 为上机编程方便，改写为 y p yn hf ( xn, yn) yq yn hf ( xn h, y p) yn 1 1 ( y p yq ) 2
y x y ,0 x 1 y(0) 0
5
解： 易知精确解为y x e x 1,
Euler公式的具体形式为yn+1 =yn +h(xn -yn ),(n=0,1, )
若取h=0.1,N=10,可由初值y0=0逐次算出.见下表
xn 0 0.1 0.2 0.3 0.9 1.0
f ( x, y) f ( x, y ) L y y ,(10.2)
由常微分方程理论知，初值问题（10.1） 有解且唯一。
常微分方程的数值解法全文
第8章常微分方程的数值解法8.4单步法的收敛性与稳定性8.4.1相容性与收敛性上面所介绍的方法都是用离散化的方法，将微分方程初值问题化为差分方程初值问题求解的.这些转化是否合理？即当h →∞时，差分方程是否能无限逼近微分方程，差分方程的解n y 是否能无限逼近微分方程初值问题的准确解()n y x ，这就是相容性与收敛性问题.用单步法(8.3.14)求解初值问题(8.1.1)，即用差分方程初值问题100(,,)()n n n n y y h x y h y x y ϕ+=+⎧⎨=⎩(8.4.1)的解作为问题(8.1.1)的近似解，如果近似是合理的，则应有()()(,(),)0 (0)y x h y x x y x h h hϕ+--→→(8.4.2)其中()y x 为问题(8.1.1)的精确解.因为0()()lim ()(,)h y x h y x y x f x y h→+-'==故由(8.4.2)得lim (,,)(,)h x y h f x y ϕ→=如果增量函数(,(),)x y x h ϕ关于h 连续，则有(,,0)(,)x y f x y ϕ=(8.4.3)定义8.3如果单步法的增量函数(,,)x y h ϕ满足条件(8.4.3)，则称单步法(8.3.14)与初值问题(8.1.1)相容.通常称(8.4.3)为单步法的相容条件.满足相容条件(8.4.3)是可以用单步法求解初值问题(8.1.1)的必要条件.容易验证欧拉法和改进欧拉法均满足相容性条件.一般地，如果单步法有p 阶精度（1p ≥），则其局部截断误差为[]1()()(,(),)()p y x h y x h x y x h O h ϕ++-+=上式两端同除以h ，得()()(,,)()p y x h y x x y h O h hϕ+--=令0h →，如果(,(),)x y x h ϕ连续，则有()(,,0)0y x x y ϕ'-=所以1p ≥的单步法均与问题(8.1.1)相容.由此即得各阶龙格－库塔法与初值问题(8.1.1)相容.定义8.4一种数值方法称为是收敛的，如果对于任意初值0y 及任意固定的(,]x a b ∈，都有lim () ()n h y y x x a nh →==+其中()y x 为初值问题(8.1.1)的精确解.如果我们取消局部化假定，使用某单步法公式，从0x 出发，一步一步地推算到1n x +处的近似值1n y +.若不计各步的舍入误差，而每一步都有局部截断误差，这些局部截断误差的积累就是整体截断误差.定义8.5称111()n n n e y x y +++=-为某数值方法的整体截断误差.其中()y x 为初值问题(8.1.1)的精确解，1n y +为不计舍入误差时用某数值方法从0x 开始，逐步得到的在1n x +处的近似值（不考虑舍入误差的情况下，局部截断误差的积累）.定理8.1设单步法(8.3.14)具有p 阶精度，其增量函数(,,)x y h ϕ关于y 满足利普希茨条件，问题(8.1.1)的初值是精确的，即00()y x y =，则单步法的整体截断误差为111()()p n n n e y x y O h +++=-=证明由已知，(,,)x y h ϕ关于y 满足利普希茨条件，故存在0L >，使得对任意的12,y y 及[,]x a b ∈，00h h <≤，都有1212(,,)(,,)x y h x y h L y y ϕϕ-≤-记1()(,(),)n n n n y y x h x y x h ϕ+=+，因为单步法具有p 阶精度，故存在0M >，使得1111()p n n n R y x y Mh ++++=-≤从而有111111111()()()(,(),)(,,)()(,(),)(,,)n n n n n n n p n n n n n n p n n n n n n e y x y y x y y y Mh y x h x y x h y h x y h Mh y x y h x y x h x y h ϕϕϕϕ+++++++++=-≤-+-≤++--≤+-+-1(1)p nMh hL e +≤++反复递推得11111101110(1)(1)1(1)(1)(1)(1)1(1)p p n n n p n n p n e Mh hL Mh hL e hL hL Mh hL e hL Mh hL e hL+++-+++++⎡⎤≤++++⎣⎦⎡⎤≤+++++++⎣⎦+-≤++因为00()y x y =，即00e =，又(1)n h b a +≤-，于是ln(1)1()(1)(1)b a b a hL n L b a h h hL hL e e --++-+≤+=≤所以()11()p L b a p n M e h e O h L -+⎡⎤≤-=⎣⎦推论设单步法具有p (1p ≥)阶精度，增量函数(,,)x y h ϕ在区域G ：, , 0a x b y h h ≤≤-∞<<+∞≤≤上连续，且关于y 满足利普希茨条件，则单步法是收敛的.当(,)f x y 在区域:,D a x b y ≤≤-∞<<+∞上连续，且关于y 满足利普希茨条件时，改进欧拉法，各阶龙格－库塔法的增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，因而它们都是收敛的.关于单步法收敛的一般结果是：定理8.2设增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，则单步法收敛的充分必要条件是相容性条件(8.4.3).8.4.2稳定性稳定性与收敛性是两个不同的概念，收敛性是在假定每一步计算都准确的前提下，讨论当步长0h →时，方法的整体截断误差是否趋于零的问题.而稳定性则是讨论舍入误差的积累能否对计算结果有严重影响的问题.定义8.6若一种数值方法在节点值n y 上有一个大小为δ的扰动，于以后各节点()m y m n >上产生的偏差均不超过δ，则称该方法是稳定的.我们以欧拉法为例进行讨论.假设由于舍入误差，实际得到的不是n y 而是n n n y y δ=+，其中n δ是误差.由此再计算一步，得到1(,)n n n n y y hf x y +=+把它与不考虑舍入误差的欧拉公式相减，并记111n n n y y δ+++=-，就有[]1(,)(,)1(,)n n n n n n y n nh f x y f x y hf x δδηδ+⎡⎤=+-=+⎣⎦其中y f f y∂=∂.如果满足条件1(,)1y n hf x η+≤，(8.4.4)则从n y 到1n y +的计算，误差是不增的，可以认为计算是稳定的.如果条件(8.4.4)不满足，则每步误差将增大.当0y f >时，显然条件(8.4.4)不可能满足，我们认为问题本身具有先天的不稳定性.当0y f <时，为了满足稳定性要求(8.4.4)，有时h 要很小.一般的，稳定性与方法有关，也与步长h 的大小有关，当然也与方程中的(,)f x y 有关.为简单起见，通常只考虑数值方法用于求解模型方程的稳定性，模型方程为y y λ'=(8.4.5)其中λ为复数.一般的方程可以通过局部线性化转化为模型方程，例如在(,)x y 的邻域(,)(,)(,)()(,)()x y y f x y f x y f x y x x f x y y y '==+-+-+略去高阶项，再作变量替换就得到u u λ'=的形式.对于模型方程(8.4.5)，若Re 0λ>，类似以上分析，可以认为方程是不稳定的.所以我们只考虑Re 0λ<的情形，这时不同的数值方法可能是数值稳定的或者是数值不稳定的.当一个单步法用于试验方程y y λ'=，从n y 计算一步得到1()n n y E h y λ+=(8.4.6)其中()E h λ依赖于所选的方法.因为通过点(,)n n x y 试验方程的解曲线（它满足,()n n y y y x y λ'==）为[]exp ()n n y y x x λ=-，而一个p 阶单步法的局部截断误差在()n n y x y =时有1111()()p n n n T y x y O h ++++=-=，所以有1exp()()()p n n y h E h y O h λλ+-=(8.4.7)这样可以看出()E h λ是h e λ的一个近似值.由(8.4.6)可以看到，若n y 计算中有误差ε，则计算1n y +时将产生误差()E h λε，所以有下面定义.定义8.7如果(8.4.6)式中，()1E h λ<，则称单步法(8.3.14)是绝对稳定的.在复平面上复变量h λ满足()1E h λ<的区域，称为方法(8.3.14)的绝对稳定区域，它与实轴的交称为绝对稳定区间.在上述定义中，规定严格不等式成立，是为了和线性多步法的绝对稳定性定义一致.事实上，()1E h λ=时也可以认为误差不增长.(1)欧拉法的稳定性欧拉法用于模型方程(8.4.5)，得1(1)n n y h y λ+=+，所以有()1E h h λλ=+.所以绝对稳定条件是11h λ+<，它的绝对稳定区域是h λ复平面上以(1,0)-为中心的单位圆，见图8.3.而λ为实数时，绝对稳定区间是(2,0)-.Im()h λRe()h λ2-1-O 图8.3欧拉法的绝对稳定区域(2)梯形公式的稳定性对模型方程，梯形公式的具体表达式为11()2n n n n h y y y y λλ++=++，即11212n nh y y h λλ++=-，所以梯形公式的绝对稳定区域为12112h h λλ+<-.化简得Re()0h λ<，因此梯形公式的绝对稳定区域为h λ平面的左半平面，见图8.4.特别地，当λ为负实数时，对任意的0h >，梯形公式都是稳定的.Im()h λRe()h λO 图8.4梯形公式的绝对稳定区域(3)龙格－库塔法的稳定性与前面的讨论相仿，将龙格－库塔法用于模型方程(8.4.5)，可得二、三、四阶龙格－库塔法的绝对稳定区域分别为211()12h h λλ++<23111()()126h h h λλλ+++<2341111()()()12624h h h h λλλλ++++<当λ为实数时，二、三、四阶显式龙格－库塔法的绝对稳定区域分别为20h λ-<<、2.510h λ-<<、 2.780h λ-<<.例8.5设有初值问题21010101(0)0xy y x x y ⎧'=-≤≤⎪+⎨⎪=⎩用四阶经典龙格－库塔公式求解时，从绝对稳定性考虑，对步长h 有何限制？解对于所给的微分方程有2100,(010)1f x x y xλ∂==-<≤≤∂+在区间[0,10]上，有201010max ||max51t x x λ<<==+由于四阶经典龙格－库塔公式的绝对稳定区间为 2.7850h λ-<<，则步长h 应满足00.557h <<.。
第六章 常微分方程解法
yn1 yn h ( xn , yn , h)
§6.1 概述
常微分方程数值解法所考虑的主要问 题有：
(1) 方法推导。即用什么样的途径来导出 递推格式； (2) 收敛性。即差分方程的解能否充分逼 近微分方程初值问题的解； (3) 误差传播。在递推过程中，每一步 都会产生截断误差和舍入误差，这个误 差是否对后续各步产生严重影响。
第六章 常微分方程的数值解法
§6.4 改进欧拉方法
(modified Euler’s method)
§6.4 改进欧拉方法
梯形方法比欧拉方法更精确，但是一 种隐式方法，求解方程计算量大。 实际计算中，迭代初始值yn+1可取欧拉 方程结果，迭代一次即可，这样的计算 公式叫改进欧拉法。
§6.4 改进欧拉方法
§6.1 概述 理论做了系统阐述。在代数数论领域，他引进了相 应的符号表示法及其计算法则，建立起被称为“李 普希兹代数”的超复数系。在微分几何方面，他自 1869年起对黎曼关于n维流形的度量结构的工作做 出进一步阐述和推广，开创了微分不变量理论的研 究，因此被认为是协变微分的奠基人之一。他的工 作后来被里奇有效地用于张量分析。
§6.1 概述
本章我们将学习一阶常微分方程的初 值问题的数值解：
dy f ( x, y ) dx y ( x0 ) y0 (1) (2)
一般情况下，方程(1)有无穷多个解， 式(2)是确定解的初始条件。
§6.1 概述
定义： 如果一元函数y(x)对一切 a x b 满足 (1) ( x, y( x)) 平面区域D
计算方法 (力学系本科生)
第六章 常微分方程 的数值解法 (Integration of ordinary differential equations)
§6.1 概述
常微分方程数值解法所考虑的主要问 题有：
(1) 方法推导。即用什么样的途径来导出 递推格式； (2) 收敛性。即差分方程的解能否充分逼 近微分方程初值问题的解； (3) 误差传播。在递推过程中，每一步 都会产生截断误差和舍入误差，这个误 差是否对后续各步产生严重影响。
第六章 常微分方程的数值解法
§6.4 改进欧拉方法
(modified Euler’s method)
§6.4 改进欧拉方法
梯形方法比欧拉方法更精确，但是一 种隐式方法，求解方程计算量大。 实际计算中，迭代初始值yn+1可取欧拉 方程结果，迭代一次即可，这样的计算 公式叫改进欧拉法。
§6.4 改进欧拉方法
§6.1 概述 理论做了系统阐述。在代数数论领域，他引进了相 应的符号表示法及其计算法则，建立起被称为“李 普希兹代数”的超复数系。在微分几何方面，他自 1869年起对黎曼关于n维流形的度量结构的工作做 出进一步阐述和推广，开创了微分不变量理论的研 究，因此被认为是协变微分的奠基人之一。他的工 作后来被里奇有效地用于张量分析。
§6.1 概述
本章我们将学习一阶常微分方程的初 值问题的数值解：
dy f ( x, y ) dx y ( x0 ) y0 (1) (2)
一般情况下，方程(1)有无穷多个解， 式(2)是确定解的初始条件。
§6.1 概述
定义： 如果一元函数y(x)对一切 a x b 满足 (1) ( x, y( x)) 平面区域D
计算方法 (力学系本科生)
第六章 常微分方程 的数值解法 (Integration of ordinary differential equations)
总结第6章常微分方程数值解法
偏微分方程：含有多个自变量的微分方程。
第6章 常微分方程数值解法
相关定义
阶数：微分方程中出现的未知函数的最高 阶导数的阶数。
线性微分方程：如果未知函数y及其各阶导 数都是一次的。
非线性微分方程： 一阶常微分方程
第6章 常微分方程数值解法
典型方程的解法：可分离变量法，常系数齐 次线性方程的解法，常系数非齐次线性方程 的解法 。
单步迭代法：计算yi+1时只用到xi+1, xi, yi。
代表：龙格—库塔法
多步迭代法：计算yi+1时除用到xi+1, xi, yi，
还用到xi-p, yi-p (p=1,2, …k)。 代表：亚当斯法
第6章 常微分方程数值解法
§2 欧拉法和改进的欧拉法 一、欧拉法基本思想
dy dx
f
(x, y)
y(x0 ) y0
离散结点a=x0＜x1＜x2＜…＜xn=b， h= xi-xi-1 xi=x0+ih
积分曲线上每一点(x,y)的切线的斜率y'(x)等 于函数f(x,y)在该点的值
第6章 常微分方程数值解法
y-y0=f(x0,y0)(x-x0)⇒ y1=y0+f(x0,y0)(x1-x0)=y0+f(x0,y0)h y-y1=f(x1,y1)(x-x1)⇒ y2=y1+f(x1,y1)(x2-x1)=y1+f(x1,y1)h
但能求解的常微分方程仍然是有限的，
大多数的常微分方程是不可能给出解析解的。
如 y' x2 y2,
y' y y(0) 1
第6章 常微分方程数值解法
从实际问题中归纳出来的微分方程， 通常主要依靠数值解法来解决。也就是 求在某些点上满足一定精度的近似解。
第6章 常微分方程数值解法
相关定义
阶数：微分方程中出现的未知函数的最高 阶导数的阶数。
线性微分方程：如果未知函数y及其各阶导 数都是一次的。
非线性微分方程： 一阶常微分方程
第6章 常微分方程数值解法
典型方程的解法：可分离变量法，常系数齐 次线性方程的解法，常系数非齐次线性方程 的解法 。
单步迭代法：计算yi+1时只用到xi+1, xi, yi。
代表：龙格—库塔法
多步迭代法：计算yi+1时除用到xi+1, xi, yi，
还用到xi-p, yi-p (p=1,2, …k)。 代表：亚当斯法
第6章 常微分方程数值解法
§2 欧拉法和改进的欧拉法 一、欧拉法基本思想
dy dx
f
(x, y)
y(x0 ) y0
离散结点a=x0＜x1＜x2＜…＜xn=b， h= xi-xi-1 xi=x0+ih
积分曲线上每一点(x,y)的切线的斜率y'(x)等 于函数f(x,y)在该点的值
第6章 常微分方程数值解法
y-y0=f(x0,y0)(x-x0)⇒ y1=y0+f(x0,y0)(x1-x0)=y0+f(x0,y0)h y-y1=f(x1,y1)(x-x1)⇒ y2=y1+f(x1,y1)(x2-x1)=y1+f(x1,y1)h
但能求解的常微分方程仍然是有限的，
大多数的常微分方程是不可能给出解析解的。
如 y' x2 y2,
y' y y(0) 1
第6章 常微分方程数值解法
从实际问题中归纳出来的微分方程， 通常主要依靠数值解法来解决。也就是 求在某些点上满足一定精度的近似解。
计算方法 第6章 常微分方程数值解
已知Euler格式 yn1 yn hf ( xn , yn )
h2 y( xn1 ) yn1 2 y''( xn )
即Euler格式具有一阶精度
如果令
y( xn1 ) y( xn1 ) 2h
y'( xn )
f ( xn , yn )
并假定 y( xn1 ) yn1, y( xn ) yn
常微分方程数值解
常微分方程的数值解法
§1 引 言 §2 欧拉方法 §3 龙格-库塔方法
2
§1 引 言
在工程和科学技术的实际问题中，常需要解常微 分方程。但常微分方程组中往往只有少数较简单和典 型的常微分方程（例如线性常系数常微分方程等）可 求出其解析解。对于变系数常微分方程的解析求解就 比较困难，而一般的非线性常微分方程就更不用说了。 在大多数情况下，常微分方程只能用近似法求解。这 种近似解法可分为两大类：一类是近似解析法，如级 数解法、逐次逼近法等；另一类则是数值解法，它给 出方程在一些离散点上的近似解。
yn
2
xn yn
令 h 0.1 将 x0 0, y0 1 代入Euler格式
步进计算结果见P106表5.1
第五章：常微分方程数值解
Euler值
y 1 2x
第五章：常微分方程数值解
Euler格式的误差分析
pn1
事实上Euler格式的每一步都存在误差，为了方便讨论y算( x)
d2x
dt 2 x(t
0
c
m )
x x
0
0 (t
t) 0
x(t ) x
0
0
5
3、常微分方程的数值解
可表示为：
xi x0 ih, i 1,2,, n
数值解法需要把连续性的问题加以离散化，从而求出离散
节点的数值解。
对常微分方程数值解法的基本出发点就是离散化。其数值 解法有个基本特点，它们都采用“步进式”，即求解过程顺着 节点排列的次序一步一步地向前推进，描述这类算法，要求给
yi , yi 1 , yi 2 ,, y0 出用已知信息
率场（又称方向场）。设 (t )是常微分方程的一个通解，则 (t )
所代表的曲线集合中任意一点的切线斜率均与该点在斜率场中 的值相同，即曲线沿斜率场所指示的方向伸展。初值问题的解
的函数曲线必开始于初值条件所定义的点(a,y(a))。
可以证明,如果函数在带形区域 R=a≤t≤b,-∞＜y＜∞}内连续 ，且关于y满足李普希兹(Lipschitz)条件，即存在常数L(它与t,y 无关)使得：
由于数值积分的矩形公式精度很低，以推断出欧拉(Euler) 格式的精度不会很高，应该很粗糙。
例1 用欧拉法解初值问题：
y y xy 2 y(0) 1
(0 x 0.6)
要求：取步长h=0.2 ,计算过程保留4位小数。
2 解: h=0.2, f ( x, y) y xy 欧拉迭代格式：
f (t, y1 ) f (t, y2 ) L y1 y2
对R内任意两个 y1 , y2 都成立,则一阶常微分方程的初值问题在 定义域a, b内存在唯一解。
7.2 常微分方程初值问题数值解法的基本思想
对常微分方程初值问题式的数值解法，就是要算出精确解 y(x)在区间a,b上的一系列离散节点：
可验证，向前欧拉格式的步长越小，由此迭代得到的截断 误差就越小，即精度就越高。但此时值得注意的是：步长合适 的选取非常重要！如果步长过小，不但会增加计算量还会导致 舍入误差的过度积累。因此，选择一个合适的步长，在截断误 差和舍入误差之间取得平衡，对减小向前欧拉格式的误差至关 重要！在这一点上与数值微分类似！
xi x0 ih, i 1,2,, n
数值解法需要把连续性的问题加以离散化，从而求出离散
节点的数值解。
对常微分方程数值解法的基本出发点就是离散化。其数值 解法有个基本特点，它们都采用“步进式”，即求解过程顺着 节点排列的次序一步一步地向前推进，描述这类算法，要求给
yi , yi 1 , yi 2 ,, y0 出用已知信息
率场（又称方向场）。设 (t )是常微分方程的一个通解，则 (t )
所代表的曲线集合中任意一点的切线斜率均与该点在斜率场中 的值相同，即曲线沿斜率场所指示的方向伸展。初值问题的解
的函数曲线必开始于初值条件所定义的点(a,y(a))。
可以证明,如果函数在带形区域 R=a≤t≤b,-∞＜y＜∞}内连续 ，且关于y满足李普希兹(Lipschitz)条件，即存在常数L(它与t,y 无关)使得：
由于数值积分的矩形公式精度很低，以推断出欧拉(Euler) 格式的精度不会很高，应该很粗糙。
例1 用欧拉法解初值问题：
y y xy 2 y(0) 1
(0 x 0.6)
要求：取步长h=0.2 ,计算过程保留4位小数。
2 解: h=0.2, f ( x, y) y xy 欧拉迭代格式：
f (t, y1 ) f (t, y2 ) L y1 y2
对R内任意两个 y1 , y2 都成立,则一阶常微分方程的初值问题在 定义域a, b内存在唯一解。
7.2 常微分方程初值问题数值解法的基本思想
对常微分方程初值问题式的数值解法，就是要算出精确解 y(x)在区间a,b上的一系列离散节点：
可验证，向前欧拉格式的步长越小，由此迭代得到的截断 误差就越小，即精度就越高。但此时值得注意的是：步长合适 的选取非常重要！如果步长过小，不但会增加计算量还会导致 舍入误差的过度积累。因此，选择一个合适的步长，在截断误 差和舍入误差之间取得平衡，对减小向前欧拉格式的误差至关 重要！在这一点上与数值微分类似！
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
dt (c x)2 (at y)2
由方程无法得到x(t), y(t)的解析解
需要用数值解法求解
0
Q(c,at)
P(x,y)
b
R(c,y )
cx
3
龙格—库塔方法的 MATLAB 实现
x(t) f (t, x), x(t0 ) x0 , x (x1, xn )T , f ( f1, fn )T
6
• clear,clf,shg • %Set the definied time • %ts=0:0.05:0.5; • %ts=0:0.1:1.6; • n=length(ts); • x0=[0 0]; • a=35;b=40;c=15; • opt1=odeset('RelTol',1e-6,'AbsTol',1e-9); • [t,x]=ode45(@jisi,ts,x0,opt1,a,b,c); • %a=35;b=40;c=15; • %opt1=odeset('RelTol',1e-6,'AbsTol',1e-9); • %[t,x]=ode45(@jisi,ts,x0,opt1,a,b,c); • %exact solution x1=c • y1=a*t; • %output t,x(t),y(t) and draw x(t),y(t) • [t,x,y1] • plot(t,x),grid,gtext('x(t)','FontSize',16), • gtext('y(t)','FontSize',16),pause
常微分方程数值解
1
实例1 海上缉私
海防某部缉私艇上的雷达发现正东方向c海里处有一艘走私船 正以速度a向正北方向行驶，缉私艇立即以最大速度b(>a)前往 拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始终 指向走私船。
• 建立任意时刻缉私艇位置及 航线的数学模型,并求解; • 求出缉私艇追上走私船的时间。
b(at y) (c x)2 (at y)2
x(0) 0, y(0) 0
设: 船速a=20 (海里/小时) 艇速b=40 (海里/小时) 距离c=15 (海里)
(x(t),y(t))
a
b
0
cx
求: 缉私艇的位置x(t),y(t) 缉私艇的航线y(x)
jisi.m, seajisi.m
5
%Creat the function for jisi %Let x(1)=x, x(2)=y function dx=jisi(t,x,a,b,c) s=sqrt((c-x(1))^2+(a*t-x(2))^2); %dx=[b*(c-x(1))/s;b*(a*t-x(2))/s]; dx=[1;1]; dx(1)=b*(c-x(1))/s; dx(2)=b*(a*t-x(2))/s;
9
实例1
模型的数值解
海上缉私(续)
t
x(t)
y(t) y1(t)
0
0
0
0
opt=odeset('RelTol',1e-6,
北
b
a
艇
c
船
2
实例1 海上缉私
建立坐标系如图: t=0 艇在(0, 0), 船在(c, 0); 船速a, 艇速b 时刻 t 艇位于P(x, y), 船到达 Q(c, at)
模型: dx bcos, dy bsin
dt
dt
dx
dt
b(c x) (c x)2 (at y)2
y
dy
b(at y)
• %draw y(x): the position of tatch jisi • plot(x(:,1),x(:,2),'r*'),grid • xlabel('x','FontSize',16),ylabel('y','FontSize',16)
7
实例1 海上缉私(续) 模型的数值解
a=20, b=40, c=15
t=0.5时缉私艇追上走私船 0.45 14.7451
0.50 15.0046
y(t) 0 0.0698 0.2924 0.6906 1.2899 2.1178 3.2005 4.5552 6.1773 8.0273 9.9979
y1(t) 0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0
t
x(t)
走私船的位置 x1(t)= c=15 y1(t)=at=20t
缉私艇的航线y(x)
10 y
0
0
0.05 1.9984
0.10 3.9854
0.15 5.9445
0.20 7.8515
8
0.25 9.6705
6
0.30 11.3496
4
2
x
0
0
5
10
15
20
0.35 12.8170 0.40 13.9806
0.4 12.5384 8.2755
0.5 13.7551 12.0830
1.2 14.9986 40.0164
1.3 14.9996 44.0165
1.4 15.0117 48.0183
1.5 15.0023 52.0146
1.6 14.9866 55.9486
y1(t) 0 3.5 7.0
10.5 14.0 17.5 42.0 45.5 49.0 52.5 56.0
10.0
8
实例1 海上缉私(续)
设b，c不变，a变大 为30, 35, 接近40,
观察解的变化 :
a=35, b=40, c=15
t=? 缉私艇追上走私船
累积误差较大 提高精度!
模型的数值解
t
x(t)
y(t)
0
0
0
0.1 3.9561 0.5058
0.2 7.5928 2.1308
0.3 10.5240 4.8283
[t,x]=ode23(@f,ts,x0,opt) 3级2阶龙格-库塔公式
[t,x]=ode45(@f,ts,x0,opt) 5级4阶龙格-库塔公式
f是待解方程写成 的函数m文件：
function dx=f(t,x) dx=[f1; f2;; fn];
ts = [t0,t1, …,tf] 输出指定时刻 t0,t1, …,tf 的函数值
ts = t0:k:tf
输出 [t0,tf] 内等分点处的函数值
x0为函数初值(n维) 输出t=ts, x为相应函数值(n维)
opt为选项，缺省时精度为：相对误差10-3， 绝对误差10-6, 计算步长按精度要求自动调整
4
实例1 海上缉私(续)
模型的数值解
ydxdt源自dydtb(c x) (c x)2 (at y)2