常微分方程的数值解
常微分方程的数值解
f ( x, y1 ) f ( x, y2 ) L y1 y2
（其中 L 为 Lipschitz 常数）则初值问题（ 1 ）存 在唯一的连续解。
求问题（1）的数值解，就是要寻找解函数在一 系列离散节点x1 < x2 <……< xn < xn+1 上的近似 值y1, y 2,…,yn 。 为了计算方便，可取 xn=x0+nh，(n=0,1,2,…)， h称为步长。
(1)，(2)式称为初值问题，(3)式称为边值问题。 在实际应用中还经常需要求解常微分方程组：
f1 ( x, y1 , y2 ) y1 ( x0 ) y10 y1 (4) f 2 ( x, y1 , y2 ) y2 ( x0 ) y20 y2
本章主要研究问题(1)的数值解法，对(2)~(4)只 作简单介绍。
得 yn1 yn hf ( xn1 , yn1 )
上式称后退的Euler方法，又称隐式Euler方法。 可用迭代法求解
二、梯形方法 由
y( xn1 ) y( xn )
xn1 xn
f ( x, y( x))dx
利用梯形求积公式： x h x f ( x, y( x))dx 2 f ( xn , y( xn )) f ( xn1 , y( xn1 ))
常微分方程的数言 简单的数值方法 Runge-Kutta方法 一阶常微分方程组和高阶方程
引言
在高等数学中我们见过以下常微分方程：
y f ( x, y, y) a x b y f ( x, y ) a x b (2) (1) (1) y ( x ) y , y ( x ) y 0 0 0 0 y ( x0 ) y0 y f ( x, y, y) a x b (3) y(a) y0 , y(b) yn
数值分析常微分方程数值解法
第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
求常微分方程的数值解
求常微分方程的数值解一、背景介绍常微分方程（Ordinary Differential Equation，ODE）是描述自然界中变化的数学模型。
常微分方程的解析解往往难以求得，因此需要寻找数值解来近似地描述其行为。
求解常微分方程的数值方法主要有欧拉法、改进欧拉法、龙格-库塔法等。
二、数值方法1. 欧拉法欧拉法是最简单的求解常微分方程的数值方法之一。
它基于导数的定义，将微分方程转化为差分方程，通过迭代计算得到近似解。
欧拉法的公式如下：$$y_{n+1}=y_n+f(t_n,y_n)\Delta t$$其中，$y_n$表示第$n$个时间步长处的函数值，$f(t_n,y_n)$表示在$(t_n,y_n)$处的导数，$\Delta t$表示时间步长。
欧拉法具有易于实现和理解的优点，但精度较低。
2. 改进欧拉法（Heun方法）改进欧拉法又称Heun方法或两步龙格-库塔方法，是对欧拉法进行了精度上提升后得到的一种方法。
它利用两个斜率来近似函数值，并通过加权平均来计算下一个时间步长处的函数值。
改进欧拉法的公式如下：$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\Delta t,y_n+k_1\Delta t)$$$$y_{n+1}=y_n+\frac{1}{2}(k_1+k_2)\Delta t$$改进欧拉法比欧拉法精度更高，但计算量也更大。
3. 龙格-库塔法（RK4方法）龙格-库塔法是求解常微分方程中最常用的数值方法之一。
它通过计算多个斜率来近似函数值，并通过加权平均来计算下一个时间步长处的函数值。
RK4方法是龙格-库塔法中最常用的一种方法，其公式如下：$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_1\Delta t}{2})$$ $$k_3=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_2\Delta t}{2})$$ $$k_4=f(t_n+\Delta t,y_n+k_3\Delta t)$$$$y_{n+1}=y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)\Delta t$$三、数值求解步骤对于给定的常微分方程，可以通过以下步骤求解其数值解：1. 确定初值条件：确定$t=0$时刻的函数值$y(0)$。
常微分方程的数值解算法
常微分方程的数值解算法常微分方程的数值解算法是一种对常微分方程进行数值计算的方法，这可以帮助我们更好地理解和研究自然现象和工程问题。
在本文中，我们将介绍一些常用的数值解算法，探讨它们的优缺点和适用范围。
常微分方程（ODE）是描述自然现象和工程问题的重要数学工具。
然而，对于许多ODE解析解是无法求出的，因此我们需要通过数值方法对其进行求解。
常微分方程可以写作：y' = f(t, y)其中，y是函数，f是给定的函数，表示y随t的变化率。
这个方程可以写成初始值问题（IVP）的形式：y'(t) = f(t,y(t))，y(t0) = y0其中，y（t0）=y0是方程的初始条件。
解决IVP问题的典型方法是数值方法。
欧拉方法欧拉方法是最简单的一阶数值方法。
在欧拉方法中，我们从初始条件开始，并在t = t0到t = tn的时间内，用以下公式逐步递推求解：y n+1 = y n + hf (t n, y n)其中，f（t n，y n）是点（t n，y n）处的导数， h = tn - tn-1是时间间隔。
欧拉方法的优点是简单易懂，容易实现。
然而，它的缺点是在整个时间段上的精度不一致。
程度取决于使用的时间间隔。
改进的欧拉方法如果我们使用欧拉方法中每个时间段的中间点而不是起始点来估计下一个时间点，精度就会有所提高。
这个方法叫做改进的欧拉方法（或Heun方法）。
公式为：y n+1 = y n + h½[f(t n, y n)+f(tn+1, yn + h f (tn, yn))]这是一个二阶方法，精度比欧拉方法高，但计算量也大一些。
对于易受噪声干扰的问题，改进的欧拉方法是个很好的选择。
Runge-Kutta方法Runge-Kutta方法是ODE计算的最常用的二阶和高阶数值方法之一。
这个方法对定义域内的每个点都计算一个导数。
显式四阶Runge-Kutta方法（RK4）是最常用的Runge-Kutta方法之一，并已得到大量实践的验证。
常微分方程数值解法
介绍常微分方程数值解法常微分方程（ordinary differential equations，ODE）可用于描述许多日常存在的物理系统。
处理ODE问题常常被称为数值求解法，这指的是找到概括ODE或者其他适用于数学模型的解决方案来模括这些ODE。
这种解决方案可能在一系列不同方案中发挥重要作用，以此来提供更好的解释和预测。
常微分方程与几何图形更为相关，它利用二维或者三维空间中曲线的绘制以及分析。
通过引入一些不同的方法，可以对不同的常微分方程中的量进行描述，使得可以通过数值方法的解析来进行研究。
数值解法可能是时间消耗较多的，但有助于验证几何图形中的某些过程，以此帮助揭示数学模型。
四种常见的常微分方程数值解法四种常见的常微分方程数值解法是：前向差分法、向后差分法、中点法和全分方法。
•前向差分法：前向差分法的基本概念是利用ODE的特定解来表达时间步的影响。
这是一种基本的数值法，可以在ODE中确定任意位置的点作为终点。
在这里，任何这样的点都可以表示为ODE右边的时间步。
•向后差分法：它是反过来基于前向差分法。
它要求对ODE中的时间步进行逆向推导，以获得某一特定点的解。
向后差分法要求推导反向解中点，以便可以从每一步中获取该点的解。
•中点法：这是一种非常基本的数值解法，可以用来求解ODE中的某一步的解，但不具有直观的方法解释。
主要的思想是在每一次时间步中通过求出ODE的中点来寻找解。
•全分方法：这是一种更复杂的数值解法，它要求将ODE中的每一步解细分并解决。
与前面提到的三种解法不同，它首先要求将ODE分解成若干离散区间，然后计算每一段区间中的点。
这种解法可以更准确地进行处理，但时间消耗较多，因此比较少被使用。
优化方案在需要解决常微分方程时，为了得到最佳的结果，有必要考虑一些优化措施。
•首先，应考虑将一个复杂的ODE拆分成一些更易解决的问题。
这样做的结果是，预见到解决此ODR的总耗时将会降低。
•其次，为了加快计算速度，可以考虑使用预解算法。
常微分方程的数值解法
常微分方程的数值解法1. 引言常微分方程是自变量只有一个的微分方程，广泛应用于自然科学、工程技术和社会科学等领域。
由于常微分方程的解析解不易得到或难以求得，数值解法成为解决常微分方程问题的重要手段之一。
本文将介绍几种常用的常微分方程的数值解法。
2. 欧拉方法欧拉方法是最简单的一种数值解法，其具体步骤如下：- 将自变量的区间等分为n个子区间；- 在每个子区间上假设解函数为线性函数，即通过给定的初始条件在每个子区间上构造切线；- 使用切线的斜率（即导数）逼近每个子区间上的解函数，并将其作为下一个子区间的初始条件；- 重复上述过程直至达到所需的精度。
3. 改进的欧拉方法改进的欧拉方法是对欧拉方法的一种改进，主要思想是利用两个切线的斜率的平均值来逼近每个子区间上的解函数。
具体步骤如下： - 将自变量的区间等分为n个子区间；- 在每个子区间上构造两个切线，分别通过给定的初始条件和通过欧拉方法得到的下一个初始条件；- 取两个切线的斜率的平均值，将其作为该子区间上解函数的斜率，并计算下一个子区间的初始条件；- 重复上述过程直至达到所需的精度。
4. 二阶龙格-库塔方法二阶龙格-库塔方法是一种更为精确的数值解法，其基本思想是通过近似计算解函数在每个子区间上的平均斜率。
具体步骤如下： - 将自变量的区间等分为n个子区间；- 在每个子区间上计算解函数的斜率，并以该斜率的平均值近似表示该子区间上解函数的斜率；- 利用该斜率近似值计算下一个子区间的初始条件，并进一步逼近解函数；- 重复上述过程直至达到所需的精度。
5. 龙格-库塔法（四阶）龙格-库塔法是目前常用的数值解法之一，其精度较高。
四阶龙格-库塔法是其中较为常用的一种，其具体步骤如下：- 将自变量的区间等分为n个子区间；- 在每个子区间上进行多次迭代计算，得到该子区间上解函数的近似值；- 利用近似值计算每个子区间上的斜率，并以其加权平均值逼近解函数的斜率；- 计算下一个子区间的初始条件，并进一步逼近解函数；- 重复上述过程直至达到所需的精度。
常微分方程数值解
常微分方程数值解常微分方程数值解是数学中的一门重要学科，主要研究如何求解常微分方程，在科学计算中有着重要的应用。
常微分方程模型是自然界中广泛存在的现象描述方法，有着广泛的应用领域。
比如，在物理学中，运动中的物体的位置、速度和加速度随时间的关系就可以通过微分方程描述；在经济学中，经济变化随时间的变化也可以用微分方程来描述。
而常微分方程数值解的求解方法则提供了一种快速、高效的计算手段。
一、常微分方程数值解的基本概念常微分方程就是一个描述自变量（通常是时间）与其导数之间关系的方程。
其一般形式如下：$\frac{dy}{dt} = f(y,t)$其中 $f(y,t)$ 是一个已知的函数。
常微分方程数值解就是对于一个常微分方程，对其进行数字计算求解的方法。
常微分方程数值解常使用数值积分的方法来求解。
由于常微分方程很少有解析解，因此数值解的求解方法显得尤为重要。
二、常微分方程数值解的求解方法常微分方程数值解的求解方法很多，以下介绍其中两种方法。
1.欧拉法欧拉法是最简单的一种数值算法，其思想是通过将一个微分方程转化为一个数值积分方程来求解。
其数值积分方程为：$y_{i+1}=y_i+hf(y_i,t_i)$其中 $h$ 为步长，可以理解为每次计算的间隔。
欧拉法的主要缺点是其精度比较低，收敛速度比较慢。
因此，当需要高精度的数值解时就需要使用其他的算法。
2.级数展开方法级数展开法是通过将一个待求解的微分方程进行Taylor级数展开来求解。
通过对Taylor级数展开的前若干项进行求和，可以得到微分方程与其解的近似解。
由于级数展开法的收敛速度很快，因此可以得到相对较高精度的数值解。
但是，当级数过多时，会出现截断误差。
因此，在实际应用中需要根据所需精度和计算资源的限制来选择适当的级数。
三、常微分方程数值解的应用常微分方程数值解在现代科学技术中有着广泛的应用。
以下介绍其中两个应用领域。
1.物理建模常微分方程的物理建模是常见的应用领域。
常微分方程的数值解法
主要内容
§1、引言 §2、初值问题的数值解法--单步法 §3、龙格-库塔方法 §4、收敛性与稳定性 §5、初值问题的数值解法―多步法 §6、方程组和刚性方程 §7、习题和总结
§1、 引 言 主要内容 ➢研究的问题 ➢数值解法的意义
1.什么是微分方程 ? 现实世界中大多数事物
使得对任意的x [a,b]及y1, y2都成立
则称 f (x,y) 对y 满足李普希兹条件，L 称为 Lipschitz常数.
就可保证方程解的存在唯一性
若 f (x,y) 在区域 G连续，关于y
满足李普希兹 条件
一阶常微分方程的初值问题的解存在且唯一. 我们以下的讨论，都在满足上述条件下进行.
一阶常微分方程组常表述为：
y(x0
)
y0
(1.2)
种 数 值 解
法
其中f (x,y)是已知函数，(1.2)是定解条件也称为 初值条件。
常微分方程的理论指出：
当 f (x,y) 定义在区域 G=(a≤x≤b,｜y｜＜∞)
若存在正的常数 L 使：
(Lipschitz)条件
| f (x, y1) f (x, y2) | L | y1 y2 | (1.3)
节点 xi a ihi,一般取hi h( (b a) / n)即等距
要计算出解函数 y(x) 在一系列节点
a x0 x1 xn b
处的近似值 yi y(xi )
y f (x, y)
y
(
x0
)
y0
a xb
(1.1) (1.2)
对微分方程(1.1)两端从 xn到xn1 进行积分
内部联系非常复杂
其状态随着 时间、地点、条件 的不同而不同
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验4 常微分方程的数值解【实验目的】1．掌握用MATLAB软件求微分方程初值问题数值解的方法；2．通过实例用微分方程模型解决简化的实际问题；3．了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。
【实验内容】题3小型火箭初始重量为1400kg，其中包括1080kg燃料。
火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。
设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
模型及其求解火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=9.8m/s2。
在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。
因此有如下二式：a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8dh/dt=v又知初始时刻t=0，v=0，h=0。
记x(1)=h，x(2)=v，根据MATLAB 可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。
程序如下：function [ dx ] = rocket( t,x )a=[(32000-0.4*x(2)^2)/(1400-18*t)]-9.8;dx=[x(2);a];endts=0:1:60;x0=[0,0];[t,x]=ode45(@rocket,ts,x0);h=x(:,1);v=x(:,2);a=[(32000-0.4*(v.^2))./(1400-18*t)]-9.8; [t,h,v,a];数据如下：t h v a0 0 0 13.061.00 6.57 13.19 13.302.00 26.44 26.58 13.453.00 59.76 40.06 13.504.00 106.57 53.54 13.435.00 166.79 66.89 13.266.00 240.27 80.02 12.997.00 326.72 92.83 12.618.00 425.79 105.22 12.159.00 536.99 117.11 11.6210.00 659.80 128.43 11.0211.00 793.63 139.14 10.3812.00 937.85 149.18 9.7113.00 1091.79 158.55 9.0214.00 1254.71 167.23 8.3315.00 1425.93 175.22 7.6516.00 1604.83 182.55 6.9917.00 1790.78 189.22 6.3618.00 1983.13 195.27 5.7619.00 2181.24 200.75 5.2120.00 2384.47 205.70 4.6921.00 2592.36 210.18 4.2222.00 2804.52 214.19 3.7923.00 3020.56 217.79 3.4124.00 3240.08 221.01 3.0725.00 3462.65 223.92 2.7726.00 3687.88 226.56 2.5027.00 3915.58 228.97 2.2728.00 4145.60 231.14 2.0629.00 4377.76 233.11 1.8930.00 4611.86 234.91 1.7431.00 4847.68 236.57 1.6232.00 5085.02 238.14 1.5133.00 5323.85 239.61 1.4134.00 5564.11 240.99 1.3335.00 5805.77 242.28 1.2736.00 6048.72 243.50 1.2137.00 6292.87 244.68 1.1738.00 6538.11 245.83 1.1339.00 6784.48 246.96 1.0940.00 7031.96 248.05 1.0741.00 7280.54 249.10 1.0542.00 7530.19 250.12 1.0343.00 7780.85 251.14 1.0244.00 8032.49 252.15 1.0045.00 8285.12 253.16 0.9946.00 8538.75 254.15 0.9847.00 8793.39 255.12 0.9748.00 9049.01 256.07 0.9749.00 9305.58 257.03 0.9650.00 9563.08 257.99 0.9551.00 9821.52 258.95 0.9452.00 10080.93 259.90 0.9353.00 10341.30 260.83 0.9354.00 10602.62 261.75 0.9455.00 10864.86 262.67 0.9456.00 11127.98 263.61 0.9357.00 11392.04 264.54 0.9158.00 11657.03 265.46 0.9159.00 11922.96 266.35 0.9260.00 12189.78 267.26 0.92因此，在引擎关闭的瞬间，火箭的速度为267.26m/s，高度为12189.78m，加速度为0.92m/s2。
（2）在第二个阶段，火箭的重量保持不变，没有向上的推力，只收到重力和空气阻力。
因此有如下关系式：a=dv/dt=(-mg-0.4v2)/m=-0.4v2/320-9.8dh/dt=v假设在80秒之前达到最高点，以60秒时的速度、高度、加速度为初始值进行计算，程序如下：function [ dx ] = rocket2( t,x )dx=[x(2);-0.4*x(2)^2/320-9.8];endts2=60:1:80;x1=[12189.78,267.26];[t2,x2]=ode45(@rocket2,ts2,x1);h2=x2(:,1);v2=x2(:,2);a2=-0.4*v2.^2./320-9.8;[t2,h2,v2,a2];数据如下：t2 h2 v2 a260.00 12189.78 267.26 -99.0861.00 12416.32 192.70 -56.2262.00 12584.73 147.43 -36.9763.00 12715.67 116.11 -26.6564.00 12819.59 92.79 -20.5665.00 12902.81 74.29 -16.7066.00 12969.22 58.96 -14.1567.00 13021.42 45.73 -12.4168.00 13061.17 33.95 -11.2469.00 13089.63 23.12 -10.4770.00 13107.61 12.91 -10.0171.00 13115.55 3.02 -9.8172.00 13113.66 -6.80 -9.8673.00 13101.90 -16.78 -10.1574.00 13079.96 -27.19 -10.7275.00 13047.26 -38.34 -11.6476.00 13002.90 -50.62 -13.0077.00 12945.47 -64.56 -15.0178.00 12872.96 -80.96 -17.9979.00 12782.31 -101.08 -22.5780.00 12668.88 -127.03 -29.97可以看到：在第60秒瞬间，加速度发生了突变，从0.92m/s2突变为-99.08m/s2；而在第71秒至第72秒之间，速度从正变为负，即速度反向，说明在第71秒中的某个时刻速度为零，火箭达到了最高点。
因此需要对这个时间段进行分析，并且找到速度减小到0的时刻和此时的高度。
以0.1为步长，在71s到72s中重新求解微分方程的数值解。
71.00 13115.16 2.98 -9.8171.10 13115.41 2.00 -9.8171.20 13115.56 1.02 -9.8071.30 13115.61 0.04 -9.8071.40 13115.57 -0.94 -9.8071.50 13115.42 -1.92 -9.80可见在t=71.3时，速度为0.04，可视为速度为零点，此时最大高度为13115.61，加速度-9.80。
综合（1），（2），可以绘出高度，速度和加速度随时间的变化曲线。
plot(t,h,t2,h2),xlabel('t/s'),ylabel('h/m'),title('高度随时间变化曲线');plot(t,v,t2,v2),xlabel('t/s'),ylabel('v/(m/s)'),title('速度随时间变化曲线');aa=[a',a2']'；tt=[t',t2']'；plot(tt,aa),xlabel('t/s'),ylabel('a/(m/s2)'),title('加速度随时间的变化曲线');题6一只小船渡过宽为d 的河流，目标是起点A 正对着的另一岸B 点。
已知河水流速v 1与船在静水中的速度v 2之比为k 。
（1）建立描述小船航线的数学模型，求其解析解；（2）设d=100m ，v 1=1m/s ，v 2=2m/s ，用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较；（3）若流速v 1=0，0.5，1.5，2m/s ，结果将如何。
模型及其求解（1）.假设在航行过程中，人们不知道水流的速度，小船的方向始终指向目标B ，因此若以B 为原点，我们可以得到如下方程组：解初值为（x, y ）=( 0, -d) 的常微分方程组，得到解析解为：其中 c=–1/d=–0.01,故有1v x dxv dt v y dy dt ⎧=-⎪⎪⎨⎪=-⎪⎩1111022k k kkx c ycy+--+-=1111()(0.01)(0.01)22kk kkx y yy+--=--+-事实上，若用matlab中计算微分方程的语句：[x,y]=dsolve('Dx=v1-v2*x/sqrt(x^2+y^2)','Dy=-v2*y/sqrt(x^2+y^2)','x(0) =0','y(0)=-d','t');则显示“Warning: Explicit solution could not be found.”即无法得到x,y关于t的分立解。