常微分方程数值方法
第九章 常微方程数值解法
第8章 序
许多科学技术问题，例如天文学中的星体运动， 许多科学技术问题，例如天文学中的星体运动，空间 技术中的物体飞行，自动控制中的系统分析， 技术中的物体飞行，自动控制中的系统分析，力学中的振 动，工程问题中的电路分析等，都可归结为常微分方程的 工程问题中的电路分析等， 初值问题。 初值问题。 所谓初值问题， 所谓初值问题，是函数及其必要的导数在积分的起始 点为已知的一类问题，一般形式为： 点为已知的一类问题，一般形式为：
⇒ y n +1 = y n −1 + 2hf ( xn , y n )
第9章 常微分方程数值解法
(8 - 4)
8-10
Euler公式的推导（ Euler公式的推导（续5） 公式的推导
上对y )=f 四、利用数值积分公式：在[xn,xn+1]上对y′(x)=f (x,y(x)) 积分 利用数值积分公式：
x0 < x1 < L < xn < L
（i=1,2,…,n）构造插值函数作为近似函数。上述离散点 i=1,2,…,n）构造插值函数作为近似函数。 相 邻两点间的距离h 称为步长， 邻两点间的距离hi=xi-1-xi 称为步长，若hi 都相等为一定数 h, 则称为定步长，否则为变步长。( x, y ( x)) 则称为定步长，否则为变步长。 a≤ x≤b y ′( x) = f 本章重点讨论如下 y (a ) = y0 一阶微分方程： 一阶微分方程： 在此基础上介绍一阶微分方程组与 8-5 第9章 常微分方程数值解法 高阶微分方程的数值解法。 高阶微分方程的数值解法。
⇒ yn +1 = yn + hf ( xn , yn ) + E ( xn , h) ⇒ yn +1 = yn + hf ( xn , yn )
常微分方程的数值解法及其应用研究
常微分方程的数值解法及其应用研究引言：常微分方程是数学中的重要分支，广泛应用于自然科学、工程技术和社会经济等领域。
常微分方程的解析解往往难以获得，因此数值解法的研究成为解决实际问题的有效手段。
本文将介绍常微分方程的数值解法以及其在各个领域的应用。
一、常微分方程的数值解法1. 欧拉方法欧拉方法是最基本的数值解法之一，通过将微分方程中的函数进行逐步的线性近似，得到方程的递推关系，并根据该关系逼近解析解。
欧拉方法具有简单、易于实现的优点，但在稳定性和精度方面存在一定的局限性。
2. 改进的欧拉方法改进的欧拉方法通过使用中点梯形公式，对欧拉方法的误差进行修正，提高了数值解的准确性。
改进的欧拉方法在简单性和准确性方面取得了一定的平衡。
3. 4阶龙格-库塔法4阶龙格-库塔法是一类常用的数值解法，通过计算多个近似解，并按照一定的权重进行加权平均，得到更高精度的数值解。
4阶龙格-库塔法具有高精度和较好的稳定性，被广泛应用于各个领域。
4. 多步法多步法是一类基于历史步长的数值解法，利用之前计算的步长来估计下一个步长的近似值。
常见的多步法包括亚当斯方法和预报校正方法等。
多步法在一定程度上提高了数值解的稳定性和准确性。
5. 常微分方程的辛方法辛方法是一类特殊的数值解法，能够保持微分方程的守恒性质。
辛方法在长时间积分和保持能量守恒方面具有优势，被广泛应用于天体力学和分子动力学等领域。
二、常微分方程数值解法的应用1. 物理科学中的应用常微分方程的数值解法在物理学中有广泛的应用，如天体力学中的行星轨道计算、量子力学中的薛定谔方程求解等。
数值解法处理了复杂的物理现象，为物理学研究提供了可行的途径。
2. 工程技术中的应用常微分方程的数值解法在工程技术中被广泛应用，如电路分析、结构力学、流体力学等。
通过数值解法，可以模拟和分析复杂的工程问题，提供设计和优化方案。
3. 经济学中的应用经济学中的许多问题可以转化为常微分方程的形式，如经济增长模型、市场供需关系等。
求常微分方程的数值解
求常微分方程的数值解一、背景介绍常微分方程（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)$。
常微分方程中的数值方法
常微分方程中的数值方法常微分方程是数学中的一个重要分支。
它主要研究的对象是随时间变化的函数。
在实际应用中，我们需要求解这些函数的解析解，但通常情况下，解析解并不容易得到，甚至是不可能得到。
因此，我们需要使用数值方法来求解这些函数的数值近似解。
在本文中，我们将介绍常微分方程中的数值方法。
一、欧拉法欧拉法是常微分方程数值解法中最基本的一种方法。
它是根据欧拉公式推导而来的。
具体地，我们可以将一阶常微分方程dy/dt=f(t,y)写成如下形式：y(t+h)=y(t)+hf(t,y(t))其中，h是步长，f(t,y)是t时刻y的导数。
欧拉法就是通过上面的公式进行逐步逼近，然后得到最终的数值解。
欧拉法的计算过程非常简单，但所得到的解可能会出现误差。
这是因为欧拉法忽略了f(t+h,y(t+h))和f(t,y(t))之间的变化。
因此，我们需要使用更为精确的数值方法来解决这个问题。
二、改进欧拉法为了解决欧拉法中的误差问题，我们可以使用改进欧拉法。
改进欧拉法又称作四阶龙格-库塔法。
它的基本思想是对欧拉法公式进行改进，以提高计算精度。
具体地，根据龙格-库塔公式，可将改进欧拉法表示为：y(t+h)=y(t)+1/6(k1+2k2+2k3+k4)其中，k1=h*f(t,y)k2=h*f(t+h/2,y+k1/2)k3=h*f(t+h/2,y+k2/2)k4=h*f(t+h,y+k3)改进欧拉法的计算过程比欧拉法要复杂些，但所得到的数值解比欧拉法更精确。
这种方法适用于一些特殊的问题，但在求解一些更为复杂的问题时，还需要使用其他的数值方法。
三、龙格-库塔法龙格-库塔法是求解常微分方程中数值解的常用方法之一。
它最常用的是四阶龙格-库塔法。
这种方法的基本思想是使用四个不同的斜率来计算数值解。
具体地，我们可以将四阶龙格-库塔法表示为：y(t+h)=y(t)+1/6(k1+2k2+2k3+k4)其中，k1=h*f(t,y)k2=h*f(t+h/2,y+k1/2)k3=h*f(t+h/2,y+k2/2)k4=h*f(t+h,y+k3)与改进欧拉法相比，龙格-库塔法的计算复杂度更高，但所得到的数值解更为精确。
常微分方程的数值解法
常微分方程的数值解法1. 引言常微分方程是自变量只有一个的微分方程，广泛应用于自然科学、工程技术和社会科学等领域。
由于常微分方程的解析解不易得到或难以求得，数值解法成为解决常微分方程问题的重要手段之一。
本文将介绍几种常用的常微分方程的数值解法。
2. 欧拉方法欧拉方法是最简单的一种数值解法，其具体步骤如下：- 将自变量的区间等分为n个子区间；- 在每个子区间上假设解函数为线性函数，即通过给定的初始条件在每个子区间上构造切线；- 使用切线的斜率（即导数）逼近每个子区间上的解函数，并将其作为下一个子区间的初始条件；- 重复上述过程直至达到所需的精度。
3. 改进的欧拉方法改进的欧拉方法是对欧拉方法的一种改进，主要思想是利用两个切线的斜率的平均值来逼近每个子区间上的解函数。
具体步骤如下： - 将自变量的区间等分为n个子区间；- 在每个子区间上构造两个切线，分别通过给定的初始条件和通过欧拉方法得到的下一个初始条件；- 取两个切线的斜率的平均值，将其作为该子区间上解函数的斜率，并计算下一个子区间的初始条件；- 重复上述过程直至达到所需的精度。
4. 二阶龙格-库塔方法二阶龙格-库塔方法是一种更为精确的数值解法，其基本思想是通过近似计算解函数在每个子区间上的平均斜率。
具体步骤如下： - 将自变量的区间等分为n个子区间；- 在每个子区间上计算解函数的斜率，并以该斜率的平均值近似表示该子区间上解函数的斜率；- 利用该斜率近似值计算下一个子区间的初始条件，并进一步逼近解函数；- 重复上述过程直至达到所需的精度。
5. 龙格-库塔法（四阶）龙格-库塔法是目前常用的数值解法之一，其精度较高。
四阶龙格-库塔法是其中较为常用的一种，其具体步骤如下：- 将自变量的区间等分为n个子区间；- 在每个子区间上进行多次迭代计算，得到该子区间上解函数的近似值；- 利用近似值计算每个子区间上的斜率，并以其加权平均值逼近解函数的斜率；- 计算下一个子区间的初始条件，并进一步逼近解函数；- 重复上述过程直至达到所需的精度。
常微分方程数值解法
第七章 常微分方程数值解法常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。
另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。
因此，有必要探讨常微分方程初值问题的数值解法。
本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。
第一节 欧拉法求解常微分方程初值问题⎪⎩⎪⎨⎧==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)公式。
以上方法称 为欧 拉法或欧拉折线法。
数值计算中的常微分方程数值模拟
数值计算中的常微分方程数值模拟在数值计算中，常微分方程（Ordinary Differential Equations，简称ODE）是一个重要的研究对象。
常微分方程的数值模拟是通过数值方法对其进行近似求解的过程，该过程对于模拟物理系统、生物学过程以及工程问题等具有重要意义。
本文将介绍常微分方程数值模拟的几种常用方法，并分析其特点与应用。
一、欧拉法（Euler's Method）欧拉法是最简单的常微分方程数值模拟方法之一，其基本思想是将连续的微分方程进行离散化，使用一阶差分近似代替微分。
具体步骤如下：1. 建立微分方程：设待求解的微分方程为dy/dx = f(x, y)，其中f(x, y)为已知函数。
2. 初始化：选择初始条件y0 = y(x0)，以及离散步长h。
3. 迭代求解：根据欧拉法的迭代公式yn+1 = yn + h * f(xn, yn)进行近似求解。
欧拉法的优点是简单易实现，但在处理复杂问题和大步长时存在精度较低的问题。
二、改进的欧拉法（Improved Euler's Method）为了提高欧拉法的精度，改进的欧拉法在迭代过程中使用两个不同的斜率近似值，从而对解进行更准确的预测并修正。
具体步骤如下：1. 建立微分方程：同欧拉法。
2. 初始化：同欧拉法。
3. 迭代求解：根据改进的欧拉法的迭代公式yn+1 = yn + h * (k1 +k2)/2进行近似求解，其中k1 = f(xn, yn)，k2 = f(xn + h, yn + h * k1)。
改进的欧拉法在精度上优于欧拉法，但仍然不适用于高精度要求的问题。
三、龙格-库塔法（Runge-Kutta Methods）龙格-库塔法是一类常微分方程数值模拟方法，通过计算多个不同次数的斜率来逼近解。
其中，四阶龙格-库塔方法是最常用的一种方法。
具体步骤如下：1. 建立微分方程：同欧拉法。
2. 初始化：同欧拉法。
3. 迭代求解：根据四阶龙格-库塔方法的迭代公式yn+1 = yn + h * (k1 + 2k2 + 2k3 + k4)/6进行近似求解，其中k1 = f(xn, yn)，k2 = f(xn + h/2, yn + h/2 * k1)，k3 = f(xn + h/2, yn + h/2 * k2)，k4 = f(xn + h, yn + h * k3)。
常微分方程数值解法
ρ ρ
n+1 n
≤1
三、梯形公式
由 分 径 y ( xn+1) = y ( xn) + 积 途 ： xn+1
∫
f ( x, y)dt
（
积分 梯形 式 且令：yn+1 = y( xn+1), yn = y( xn) 用 公 ， h 则 yn+1 = yn + ( f (xn , yn) + f (xn+1 , yn+1)) 得： 2
第九章 常微分方程数值解法
§1 、引言
一 常 分 程 初 问 : 阶 微 方 的 值 题 dy dx = f (x, y) y( x0) = y0
'
a ≤ x ≤b
2 y 例 : 方 程 xy -2 y = 4 x ⇒ y = + 4 x 2 y 令 ：f ( x , y ) = + 4 且 给 出 初 值 y (1 )= -3 x 就 得 到 一 阶 常 微 分 方 程 的 初 值 问 题 ： 2 y dy = f (x, y) = + 4 dx x y(1) = − 3
n n n n n 2 // n n+1
~
y
n+1
= yn + hf ( xn, yn ) = y(xn) + hf
n+1
~
y
n+1
( x , y( x ))
n n
则 T = y( x ) − = h y (ξ ) x y 2 ~
// n+1 n+1
2
n
< ξ < xn+1
令
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
常微分方程数值方法1、欧拉方法：1,,1,0),,(1-=+=+n k y t hf y y k k k k .function E=euler(f,a,b,ya,n)% Input - f is the function entered as a string 'f'% - a and b are the left and right end points% - ya is the initial condition y(a)% - n is the number of steps% Output - E=[T' Y'] where T is the vector of abscissas and% Y is the vector of ordinatesh=(b-a)/n;T=zeros(1,n+1);Y=zeros(1,n+1);T=a:h:b;Y(1)=ya;for j=1:nY(j+1)=Y(j)+h*feval(f,T(j),Y(j));endE=[T' Y'];【例】 用欧拉方法求解区间]3,0[内的初值问题：1)0(,2'=-=y yt y 。
f=inline('(t-y)/2','t','y');a=0;b=3;ya=1;n=12;%n=3,6,12,24,48,96...E=euler(f,a,b,ya,n),plot(E(:,1),E(:,2),'r*'),hold on符号解：y=dsolve('Dy=(t-y)/2','y(0)=1')h=(3-0)/12;t=0:h:3;y=eval(y);[t' y']用图比较数值解：(f 为ode 函数文件)ode45('f',[0,3],1)2、休恩(Huen)方法(即改进Euler 方法)：1,,1,0)],,(,(),([211-=+++=++n k y t hf y t f y t f hy y k k k k k k k kfunction H=heun(f,a,b,ya,n)% Input - f is the function entered as a string 'f'% - a and b are the left and right end points% - ya is the initial condition y(a)% - n is the number of steps% Output - H=[T' Y'] where T is the vector of abscissas and% Y is the vector of ordinatesh=(b-a)/n;T=zeros(1,n+1);Y=zeros(1,n+1);T=a:h:b;Y(1)=ya;for j=1:nk1=feval(f,T(j),Y(j));k2=feval(f,T(j+1),Y(j)+h*k1);Y(j+1)=Y(j)+(h/2)*(k1+k2);endH=[T' Y'];【例】 用休恩方法求解区间]3,0[内的初值问题：1)0(,2'=-=y yt y 。
f=inline('(t-y)/2','t','y');a=0;b=3;ya=1;n=12;H=heun(f,a,b,ya,n),plot(H(:,1),H(:,2),'r*'),hold on用图比较数值解：(f 为ode 函数文件)ode45('f',[0,3],1)3、四阶龙格－库塔方法：)22(643211k k k k hy y k k ++++=+function R=rk4(f,a,b,ya,n)% Input - f is the function entered as a string 'f'% - a and b are the left and right end points% - ya is the initial condition y(a)% - n is the number of steps% Output - R=[T' Y'] where T is the vector of abscissas and% Y is the vector of ordinatesh=(b-a)/n;T=zeros(1,n+1);Y=zeros(1,n+1);T=a:h:b;Y(1)=ya;for j=1:nk1=h*feval(f,T(j),Y(j));k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);k4=h*feval(f,T(j)+h,Y(j)+k3);Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;endR=[T' Y'];【例】 用4阶龙格－库塔方法求解区间]3,0[内的初值问题：1)0(,2'=-=y yt y 。
f=inline('(t-y)/2','t','y');a=0;b=3;ya=1;n=3;R=rk4(f,a,b,ya,n)plot(R(:,1),R(:,2),'r*'),hold on用图比较数值解：(f 为ode 函数文件)ode45('f',[0,3],1)4、阿当姆斯－巴什弗斯－摩尔顿方法：（预测－校正公式）⎪⎪⎩⎪⎪⎨⎧++-+=+-+-+=+--+---)),(9195(24)5559379(241121123p x f f f f h y y ff f f h y p k k k k k k k k k k kfunction A=adams(f,T,Y)% Input - f is the function entered as a string 'f'% - T is the vector of abscissas% - Y is the vector of ordinates% Remark:The first four coordinates of T and% Y must have starting values obtained with RK4% Output - A=[T' Y'] where T is the vector of abscissas and% Y is the vector of ordinatesn=length(T);if n<5,break,end;F=zeros(1,4);F=feval(f,T(1:4),Y(1:4));h=T(2)-T(1);for k=4:n-1% Predictorp=Y(k)+(h/24)*(F*[-9 37 -59 55]');T(k+1)=T(1)+h*k;F=[F(2) F(3) F(4) feval(f,T(k+1),p)];% CorrectorY(k+1)=Y(k)+(h/24)*(F*[1 -5 19 9]');F(4)=feval(f,T(k+1),Y(k+1));endA=[T' Y'];【例】 求解区间]3,0[内的初值问题：1)0(,2'=-=y y t y 。
f=inline('(t-y)/2','t','y'); T=zeros(1,25);Y=zeros(1,25);T=0:1/8:3; Y(1:4)=[1 0.94323919 0.89749071 0.86208736];A=adams(f,T,Y)5、米尔尼－辛普生（Milne-simpson ）方法：预测:)22(341231k k k k k f f f h y p +-+=---+ 修正：292811k k k k p y p m -+=++，),(111+++=k k k m t f f 校正:)4(31111+--++++=k k k k k f f f h y y function M=milne(f,T,Y)% Input - f is the function entered as a string 'f'% - T is the vector of abscissas% - Y is the vector of ordinates% Remark:The first four coordinates of T and% Y must have starting values obtained with RK4% Output - M=[T' Y'] where T is the vector of abscissas and% Y is the vector of ordinatesn=length(T);if n<5,break,end;F=zeros(1,4);F=feval(f,T(1:4),Y(1:4));h=T(2)-T(1);pold=0;yold=0;for k=4:n-1% Predictorpnew=Y(k-3)+(4*h/3)*(F(2:4)*[2 -1 2]');% Modifierpmod=pnew+28*(yold-pold)/29;T(k+1)=T(1)+h*k;F=[F(2) F(3) F(4) feval(f,T(k+1),pmod)];% CorrectorY(k+1)=Y(k-1)+(h/3)*(F(2:4)*[1 4 1]');pold=pnew;yold=Y(k+1);F(4)=feval(f,T(k+1),Y(k+1));endM=[T' Y'];【例】 求解区间]3,0[内的初值问题：1)0(,2'=-=y y t y 。
f=inline('(t-y)/2','t','y'); T=zeros(1,25);Y=zeros(1,25);T=0:1/8:3; Y(1:4)=[1 0.94323919 0.89749071 0.86208736];M=milne(f,T,Y)6、哈明(Hamming)方法：预测：)22(341231k k k k k f f f h y p +-+=---+ 校正：)),(2(83)9(8111121++--+++-++-=k k k k k k k p t f f f h y y y function H=hamming(f,T,Y)% Input - f is the function entered as a string 'f'% - T is the vector of abscissas% - Y is the vector of ordinates% Remark:The first four coordinates of T and% Y must have starting values obtained with RK4% Output - H=[T' Y'] where T is the vector of abscissas and% Y is the vector of ordinatesn=length(T);if n<5,break,end;F=zeros(1,4);F=feval(f,T(1:4),Y(1:4));h=T(2)-T(1);pold=0;cold=0;for k=4:n-1% Predictorpnew=Y(k-3)+(4*h/3)*(F(2:4)*[2 -1 2]');% Modifierpmod=pnew+112*(cold-pold)/121;T(k+1)=T(1)+h*k;F=[F(2) F(3) F(4) feval(f,T(k+1),pmod)];% Correctorcnew=(9*Y(k)-Y(k-2)+3*h*(F(2:4)*[-1 2 1]'))/8;Y(k+1)=Y(k-1)+(h/3)*(F(2:4)*[1 4 1]');pold=pnew;cold=cnew;F(4)=feval(f,T(k+1),Y(k+1));endH=[T' Y'];【例】 求解区间]3,0[内的初值问题：1)0(,2'=-=y y t y 。