常微分方程(组)的数值解法
重要：MATLAB常微分方程(组)数值解法
Matlab常微分方程求解问题分类
边值问题：
初值问题：
• 定解附加条件在自变量 的一端
• 一般形式为： y' f (x, y)
y(a)
y0
• 初值问题的数值解法一 般采用步进法，如 Runge-Kutta法
➢ 在自变量两端均给定附加 条件
y' f (x, y)
➢ 一般形式：y(a)y1, y(b)y2
1.根据常微分方程要求的求解精度与速度要求
求解初值问题：
y
'
y
2x y
y ( 0 ) 1
(0x1)
比较ode45和ode23的求解精度和速度
ode45和ode23的比较-1
function xODE clear all clc
format long
y0 = 1; [x1,y1] = ode45(@f,[0,1],y0); [x2,y2] = ode23(@f,[0,1],y0); plot(x1,y1,'k-',x2,y2,'b--') xlabel('x') ylabel('y')
rD = k(3)*C(2)-k(5)*C(4);
rE = k(4)*C(3)+k(5)*C(4);
% Mass balances dCdt = [rA; rB; rC; rD; rE];
三个串联的CSTR等温反应器（例4-3）
function IsothermCSTRs clear all clc CA0 = 1.8; % kmol/m^3 CA10 = 0.4; % kmol/m^3 CA20 = 0.2; % kmol/m^3 CA30 = 0.1; % kmol/m^3 k = 0.5; % 1/min tau = 2; stoptime = 2.9; % min [t,y] = ode45(@Equations,[0 stoptime],[CA10 CA20 CA30],[],k,CA0,tau); disp(' Results:') disp(' t CA1 CA2 CA3') disp([t,y]) plot(t,y(:,1),'k--',t,y(:,2),'b:',t,y(:,3),'r-') legend('CA_1','CA_2','CA_3') xlabel('Time (min)') ylabel('Concentration') % -----------------------------------------------------------------function dydt = Equations(t,y,k,CA0,tau) CA1 = y(1); CA2 = y(2); CA3 = y(3); dCA1dt = (CA0-CA1)/tau - k*CA1; dCA2dt = (CA1-CA2)/tau - k*CA2; dCA3dt = (CA2-CA3)/tau - k*CA3; dydt = [dCA1dt; dCA2dt; dCA3dt];
常微分方程组的解法
常微分方程组的解法常微分方程组是由多个关于未知函数及其导数的方程组成的方程组，它是数学中的重要研究对象。
常微分方程组的解法可以分为解析解法和数值解法两种。
解析解法是指通过数学方法求出常微分方程组的解析表达式。
常微分方程组的解析解法主要包括分离变量法、一阶线性方程法、变量代换法、常数变易法、特殊函数法等。
其中，分离变量法是指将常微分方程组中的各个变量分离出来，然后对每个变量分别积分，最后得到常微分方程组的解析解。
一阶线性方程法是指将常微分方程组转化为一阶线性方程，然后通过求解一阶线性方程来得到常微分方程组的解析解。
变量代换法是指通过合适的变量代换将常微分方程组转化为更简单的形式，然后通过求解简化后的方程组得到常微分方程组的解析解。
常数变易法是指将常微分方程组中的常数作为未知量，然后通过求解常数得到常微分方程组的解析解。
特殊函数法是指通过特殊函数的性质求解常微分方程组，如指数函数、三角函数等。
数值解法是指通过计算机数值计算的方法求出常微分方程组的数值解。
常微分方程组的数值解法主要包括欧拉法、龙格-库塔法、变步长法等。
其中，欧拉法是一种简单的数值解法，它的基本思想是将常微分方程组的解曲线上的点离散化为一系列点，然后通过计算机逐步求解得到常微分方程组的数值解。
龙格-库塔法是一种高阶数值解法，它通过计算机采用多个不同的计算公式来逼近常微分方程组的解曲线，从而得到更为准确的数值解。
变步长法是一种自适应数值解法，它通过计算机根据误差大小自动调整步长大小，从而得到更为准确的数值解。
常微分方程组的解法包括解析解法和数值解法两种，每种方法都有其适用的范围和优缺点。
在实际应用中，需要根据具体问题的性质和求解要求选择合适的解法来求解常微分方程组。
常微分方程的数值解
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
常微分方程组求解的小波-galerkin法
常微分方程组求解的小波-galerkin法一、引言常微分方程组在科学、工程、数学等领域中有着广泛的应用。
小波-Galerkin法是一种常用的常微分方程组的数值解法，它结合了小波分析和Galerkin方法的特点，具有较高的精度和稳定性。
本论文将详细介绍小波-Galerkin法的原理、算法和应用。
二、小波分析基础小波分析是一种时间-频率分析方法，它能够提供信号在不同时间和频率上的局部化信息。
通过选择合适的小波函数，可以有效地去除噪声，提取信号的边缘和特征。
在求解常微分方程组时，小波分析可以用于构造基函数，以减少数值解的误差。
三、Galerkin方法Galerkin方法是求解偏微分方程的一种数值方法，它基于有限元思想，通过将原始问题转化为一系列简单的问题，从而得到精确的数值解。
在求解常微分方程组时，Galerkin方法可以将微分方程转化为等价的积分形式，通过求解积分方程得到数值解。
四、小波-Galerkin法原理小波-Galerkin法将小波分析和Galerkin方法相结合，通过选择合适的小波基函数和有限元空间，将常微分方程组转化为一系列简单的代数方程组，从而得到精确的数值解。
该方法具有较高的精度和稳定性，适用于求解各种类型的常微分方程组。
五、算法实现小波-Galerkin法的实现主要包括以下步骤：1.选取合适的小波基函数和有限元空间；2.将常微分方程组转化为等价的积分形式；3.对方程组中的每一个方程应用Galerkin方法，得到一系列代数方程；4.对每个方程应用小波变换，提取时间-频率上的局部化信息；5.对方程进行数值求解。
六、应用举例通过具体例子，展示小波-Galerkin法在求解常微分方程组中的应用。
选择一组简单的常微分方程组，采用小波-Galerkin法进行数值求解，并与传统的方法进行比较，分析结果的可信度和精度。
七、结论本论文详细介绍了小波-Galerkin法的原理、算法和应用。
通过将小波分析和Galerkin方法相结合，该方法具有较高的精度和稳定性，适用于求解各种类型的常微分方程组。
分别利用欧拉法和改进欧拉法求解微分方程组的数值解
分别利用欧拉法和改进欧拉法求解微分方程组的数值解欧拉法（Euler’s Method）和改进欧拉法（Improved Euler’s Method），是求解常微分方程数值解的两种常用方法。
它们都属于一阶精度的显式迭代算法。
首先，我们来介绍一下欧拉法。
欧拉法是一种简单的数值求解算法，它基于微分方程的定义，将微分方程转化为差分方程。
考虑一个一阶常微分方程 dy/dx = f(x, y)，并给定初始条件 y(x0)= y0，我们希望求解在给定区间 [x0, xn] 上方程的数值解。
首先，我们将区间 [x0, xn] 平均分成 N 个小区间，每个小区间的长度为 h = (xn - x0) / N。
然后，我们可以使用以下的欧拉迭代公式计算数值解：y[i+1] = y[i] + h * f(x[i], y[i])其中，x[i] = x0 + i * h，y[i] 是在点 x[i] 处的数值解。
通过不断迭代上述公式，我们可以获得[x0, xn] 上微分方程的数值解。
欧拉法的优点在于简单易懂，计算速度较快。
然而，欧拉法的缺点是精度较低，误差随着步长h 的增大而增大。
为了提高精度，我们可以使用改进欧拉法。
改进欧拉法，也称为龙格–库塔算法（Runge-Kutta Method）或四阶龙格–库塔方法，是一种基于欧拉法的改进算法。
改进欧拉法使用了更多的近似取值，以减小误差。
与欧拉法类似，我们将区间 [x0, xn] 平均分成 N 个小区间，每个小区间的长度为 h = (xn - x0) / N。
然后，我们可以使用以下的公式计算数值解：k1 = h * f(x[i], y[i])k2 = h * f(x[i] + h/2, y[i] + k1/2)y[i+1] = y[i] + k2其中，k1 和 k2 是计算过程中的辅助变量。
通过不断迭代上述公式，我们可以获得 [x0, xn] 上微分方程的数值解。
改进欧拉法相对于欧拉法而言，计算精度更高。
常微分方程的数值解法及其应用研究
常微分方程的数值解法及其应用研究引言：常微分方程是数学中的重要分支，广泛应用于自然科学、工程技术和社会经济等领域。
常微分方程的解析解往往难以获得，因此数值解法的研究成为解决实际问题的有效手段。
本文将介绍常微分方程的数值解法以及其在各个领域的应用。
一、常微分方程的数值解法1. 欧拉方法欧拉方法是最基本的数值解法之一，通过将微分方程中的函数进行逐步的线性近似，得到方程的递推关系，并根据该关系逼近解析解。
欧拉方法具有简单、易于实现的优点，但在稳定性和精度方面存在一定的局限性。
2. 改进的欧拉方法改进的欧拉方法通过使用中点梯形公式，对欧拉方法的误差进行修正，提高了数值解的准确性。
改进的欧拉方法在简单性和准确性方面取得了一定的平衡。
3. 4阶龙格-库塔法4阶龙格-库塔法是一类常用的数值解法，通过计算多个近似解，并按照一定的权重进行加权平均，得到更高精度的数值解。
4阶龙格-库塔法具有高精度和较好的稳定性，被广泛应用于各个领域。
4. 多步法多步法是一类基于历史步长的数值解法，利用之前计算的步长来估计下一个步长的近似值。
常见的多步法包括亚当斯方法和预报校正方法等。
多步法在一定程度上提高了数值解的稳定性和准确性。
5. 常微分方程的辛方法辛方法是一类特殊的数值解法，能够保持微分方程的守恒性质。
辛方法在长时间积分和保持能量守恒方面具有优势，被广泛应用于天体力学和分子动力学等领域。
二、常微分方程数值解法的应用1. 物理科学中的应用常微分方程的数值解法在物理学中有广泛的应用，如天体力学中的行星轨道计算、量子力学中的薛定谔方程求解等。
数值解法处理了复杂的物理现象，为物理学研究提供了可行的途径。
2. 工程技术中的应用常微分方程的数值解法在工程技术中被广泛应用，如电路分析、结构力学、流体力学等。
通过数值解法，可以模拟和分析复杂的工程问题，提供设计和优化方案。
3. 经济学中的应用经济学中的许多问题可以转化为常微分方程的形式，如经济增长模型、市场供需关系等。
常微分方程的数值解法
常微分方程的数值解法1. 引言常微分方程是自变量只有一个的微分方程，广泛应用于自然科学、工程技术和社会科学等领域。
由于常微分方程的解析解不易得到或难以求得，数值解法成为解决常微分方程问题的重要手段之一。
本文将介绍几种常用的常微分方程的数值解法。
2. 欧拉方法欧拉方法是最简单的一种数值解法，其具体步骤如下：- 将自变量的区间等分为n个子区间；- 在每个子区间上假设解函数为线性函数，即通过给定的初始条件在每个子区间上构造切线；- 使用切线的斜率（即导数）逼近每个子区间上的解函数，并将其作为下一个子区间的初始条件；- 重复上述过程直至达到所需的精度。
3. 改进的欧拉方法改进的欧拉方法是对欧拉方法的一种改进，主要思想是利用两个切线的斜率的平均值来逼近每个子区间上的解函数。
具体步骤如下： - 将自变量的区间等分为n个子区间；- 在每个子区间上构造两个切线，分别通过给定的初始条件和通过欧拉方法得到的下一个初始条件；- 取两个切线的斜率的平均值，将其作为该子区间上解函数的斜率，并计算下一个子区间的初始条件；- 重复上述过程直至达到所需的精度。
4. 二阶龙格-库塔方法二阶龙格-库塔方法是一种更为精确的数值解法，其基本思想是通过近似计算解函数在每个子区间上的平均斜率。
具体步骤如下： - 将自变量的区间等分为n个子区间；- 在每个子区间上计算解函数的斜率，并以该斜率的平均值近似表示该子区间上解函数的斜率；- 利用该斜率近似值计算下一个子区间的初始条件，并进一步逼近解函数；- 重复上述过程直至达到所需的精度。
5. 龙格-库塔法（四阶）龙格-库塔法是目前常用的数值解法之一，其精度较高。
四阶龙格-库塔法是其中较为常用的一种，其具体步骤如下：- 将自变量的区间等分为n个子区间；- 在每个子区间上进行多次迭代计算，得到该子区间上解函数的近似值；- 利用近似值计算每个子区间上的斜率，并以其加权平均值逼近解函数的斜率；- 计算下一个子区间的初始条件，并进一步逼近解函数；- 重复上述过程直至达到所需的精度。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
初值问题：
• 定解附加条件在自变量 的一端
在自变量两端均给定附加 条件
y ' f ( x, y ) 一般形式： y (a) y1 , y (b) y2
• 一般形式为： y ' f ( x, y )
y (a) y 0
边值问题可能有解、也可 能无解，可能有唯一解、 也可能有无数解 边值问题有3种基本解法 • 迭加法 • 打靶法 • 松弛法
function f=fun(t,y) a=-exp(-t)+cos(2*pi*t)*exp(-2*t); b=cos(2*pi*t); f=[y(2);-a*y(2)^2-b*y(1)+exp(t)*b];
方程变形：令y1＝y；y2＝y’ 则原方程等价于：
y1' y 2 ' 2 y 2 ay2 by1 e t cos 2t
ode45和ode23的比较-1
function xODE clear all clc format long y0 = 1; [x1,y1] = ode45(@f,[0,1],y0); [x2,y2] = ode23(@f,[0,1],y0); plot(x1,y1,'k-',x2,y2,'b--') xlabel('x') ylabel('y') % -----------------------------------------------------------------function dydx = f(x,y) dydx = y - 2*x/y;
常微分方程初值问题解算指令比较
解算指令 ode45 ode23 算法 四五阶Runge-Kutta法 二三阶Runge-Kutta法 精度 较高 低
ode113 ode15s
ode23s ode23t ode23tb
可变阶Adams-Bashforth-Moulton法 基于数值差分的可变阶方法（BDFs， Gear）
二阶改进的Rosenbrock法 使用梯形规则 TR-BDF2（隐式Runge-Kutta法）
低～中
低 适中 低
ode解算指令的选择(1)
1.根据常微分方程要求的求解精度与速度要求
2x y' y y y ( 0) 1
求解初值问题：
(0 x 1)
比较ode45和ode23的求解精度和速度
高阶微分方程odefile的编写
求解： y"a(t )( y' ) 2 b(t ) y et cos 2t
a(t ) e t cos 2te 2t , b(t ) cos(2t )
y(0)=0，y'(0)=1，
本例的难度： 方程系数非线性 方程高阶，非标准形式
可在odefile中定义
常微分方程组与单个常微分方程求解方法相同， 只需在编写odefile时将整个方程组作为一个向量 输出。 function f=fun(x,y) dy1dx = 0.04*(1-y(1))-(1-y(2)).*y(1)+0.0001*(1-y(2)).^2; dy2dx = -1e4*dy1dx + 3000*(1-y(2)).^2; f = [dy1dx; dy2dx];
• 初值问题的数值解法一 般采用步进法，如 Runge-Kutta法
Matlab求解常微分方程初值问题方法
将待求解微分方程（组）转化为标准形式， “翻译”成Matlab可以理解的语言
编写odefile文件
选择合适的解算指令求解问题
根据求解问题的要求，设置解算指令的调用格式
Matlab求解初值问题函数
解算指令的使用方法
[T,Y]=ode45(@fun, TSPAN,Y0) 输出变量T为返回时间列向量；解矩阵Y的每一行对应于T的 一个元素，列数与求解变量数相等。
@fun为函数句柄，为根据待求解的ODE方程所编写的ode文
件（odefile）； TSPAN＝[T0 TFINAL]是微分系统y'＝F(t,y)的积分区间；Y0 为初始条件
ode解算指令的选择(2)
2.根据常微分方程组是否为刚性方程
y ' Ay b( x) y ( x0 ) y0
如果系数矩阵A的特征值连乘积小于零， 且绝对值最大和最小的特征值之比（刚性 比）很大，则称此类方程为刚性方程 刚性比：100/0.01＝10000
y1' 0.01 y1 99.99 y2 ' y2 100 y2 y (0) 2, y (0) 1 2 1
2.3 常微分方程（组）的数值解法
知识要点
常微分方程初值问题---ode45,0de23
微分方程在化工模型中的应用
•间歇反应器的计算 •活塞流反应器的计算
•全混流反应器的动态模拟
•定态一维热传导问题
•逆流壁冷式固定床反应器一维模型
•固定床反应器的分散模型
Matlab常微分方程求解问题分类
边值问题：
% -----------------------------------------------------------------function dydx = f(x,y) % Define simultaneous ODE equations dydx = [3*y(1) + 2*y(2); 4*y(1)+y(2)];
三个串联的CSTR等温反应器（例4-3）
function IsothermCSTRs clear all clc CA0 = 1.8; % kmol/m^3 CA10 = 0.4; % kmol/m^3 CA20 = 0.2; % kmol/m^3 CA30 = 0.1; % kmol/m^3 k = 0.5; % 1/min tau = 2; stoptime = 2.9; % min [t,y] = ode45(@Equations,[0 stoptime],[CA10 CA20 CA30],[],k,CA0,tau); disp(' Results:') disp(' t CA1 CA2 CA3') disp([t,y]) plot(t,y(:,1),'k--',t,y(:,2),'b:',t,y(:,3),'r-') legend('CA_1','CA_2','CA_3') xlabel('Time (min)') ylabel('Concentration') % -----------------------------------------------------------------function dydt = Equations(t,y,k,CA0,tau) CA1 = y(1); CA2 = y(2); CA3 = y(3); dCA1dt = (CA0-CA1)/tau - k*CA1; dCA2dt = (CA1-CA2)/tau - k*CA2; dCA3dt = (CA2-CA3)/tau - k*CA3; dydt = [dCA1dt; dCA2dt; dCA3dt];
指令 含义 指令 含义
ode23
普通2-3阶法解 ODE
普通4-5阶法解 ODE
odefile
ODE文件格式
ode45
odeset
ode113
解 算 ode23t ode15s ode23s ode23tb
普通变阶法解ODE
解适度刚性ODE 变阶法解刚性ODE 低阶法解刚性ODE 低阶法解刚性ODE
选 项
创建、更改ODE选项的 设置
odeget
读取ODE选项的设置
odefile
odefile是一个Matlab函数文件，一般作为整个求 解程序的一个子函数，表示ode求解问题 对于程序通用性要求不高的场合，只需将原有模 型写成标准形式，然后“翻译”成Matlab语言即 可
odefile的编写---常微分方程
刚性常微分方程组求解
function demo1 figure ode23s(@fun,[0,100],[0;1]) hold on, pause ode45(@fun,[0,100],[0;1]) %-------------------------------------------------------------------------function f=fun(x,y) dy1dx = 0.04*(1-y(1))-(1-y(2)).*y(1)+0.0001*(1-y(2)).^2; dy2dx = -1e4*dy1dx + 3000*(1-y(2)).^2; f = [dy1dx; dy2dx];
ode45和ode23的比较-2
function xODEs clear all clc y0 = [0 1]; [x1,y1] = ode45(@f,[0:0.2:1],y0); [x2,y2] = ode23(@f,[0:0.2:1],y0); disp('Results by using ode45()：') disp(' x y(1) y(2)') disp([x1 y1]); disp('Results by using ode23()：') disp(' x y(1) y(2)') disp([x2 y2]);
间歇反应器中串联－平行复杂反应系统（例4-1）
function BatchReactor clear all clc T = 224.6 + 273.15; % Reactor temperature, Kelvin R = 8.31434; % Gas constant, kJ/kmol K % Arrhenius constant, 1/s k0 = [5.78052E+10 3.92317E+12 1.64254E+4 6.264E+8]; % Activation energy, kJ/kmol Ea = [124670 150386 77954 111528]; % 初始浓度C0(i), kmol/m^3 C0 = [1 0 0 0 0]; tspan = [0 1e4]; [t,C] = ode45(@MassEquations, tspan, C0,[],k0,Ea,R,T); % 绘图 plot(t,C(:,1),'r-',t,C(:,2),'k:',t,C(:,3),'b-.',t,C(:,4),'k--'); xlabel('Time (s)'); ylabel('Concentration (kmol/m^3)'); legend('A','B','C','D') CBmax = max(C(:,2)); % CBmax: the maximum concentration of B, kmol/m^3 yBmax = CBmax/C0(1) % yBmax: the maximum yield of B index = find(C(:,2)==CBmax); t_opt = t(index) % t_opt: the optimum batch time, s % -----------------------------------------------------------------function dCdt = MassEquations(t,C,k0,Ea,R,T) % Reaction rate constants, 1/s k = k0.*exp(-Ea/(R*T)); k(5) = 2.16667E-04; % Reaction rates, kmoles/m3 s rA = -(k(1)+k(2))*C(1); rB = k(1)*C(1)-k(3)*C(2); rC = k(2)*C(1)-k(4)*C(3); rD = k(3)*C(2)-k(5)*C(4); rE = k(4)*C(3)+k(5)*C(4); % Mass balances dCdt = [rA; rB; rC; rD; rE];