常微分方程(组)的数值解法
重要：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、龙格-库塔方法 §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 进行积分
内部联系非常复杂
其状态随着 时间、地点、条件 的不同而不同
常微分方程数值解法
h
2
所以,只要令
2
f
xn
f n f yn O ( h )
3
1+2=1, 2=1/2, 2=1/2
(7.4)
若取=1,则得1=2=1/2,=1,此时公式(8.3)就是改进的 Euler公式; 若取1=0,则得2=1,==1/2,公式(8.3)为
y n 1 y n hK 2 K 1 f (xn , yn ) K f ( x 1 h, y 2 n n 2
……
……
f ( x n P h, y n h
i 1
p 1
pi
Ki)
其中i,i,ij为待定参数. 若此公式的局部截断误差为
O(h3),称此公式为p阶Runge-kutta方法,简称p阶R-K方法.
对于p=2的情形, 应有
K 1 f (xn , yn ) K f ( x h , y hK ) n n 1 2
显式表示出来,称这类差分公式为显式公式,而梯形公式中,
需要更多的计算量,但其计算稳定性较好.
§2 改进的Euler方法和Taylor展开方法
§2.1 改进的Euler方法 从数值积分的角度来看,梯形公式
y0
y n 1 y n h 2 , n 0 ,1, 2 , N 1 [ f ( x n , y n ) f ( x n 1 , y n 1 )]
2
1
分别取步长h=0.2 ,0.1 ,0.05,计算结果如下
h
h=0.2
h=0.1
h=0.05
xn 0.00 0.40 0.80 1.20 1.60 2.00 0.00 0.40 0.80 1.20 1.60 2.00 0.00 0.40 0.80 1.20 1.60 2.00
常微分方程数值解法
ρ ρ
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.欧拉法欧拉法是最基础的数值解法之一，它的思想是将时间连续化，将微分方程转化为差分方程。
对于一个一阶常微分方程y'=f(x,y)，其欧拉公式为：y_{n+1}=y_n+hf(x_n,y_n)其中h为步长，x_n和y_n为第n个时间点上x和y的取值。
2.改进欧拉法改进欧拉法是对欧拉法的改良，其公式如下：y_{n+1}=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},y_n+hf(x_n,y_n))] 3.四阶龙格-库塔方法四阶龙格-库塔方法是目前最常用的数值解法之一。
其公式如下：k_1=f(x_n,y_n)k_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}k_1)k_3=f(x_n+\frac{h}{2},y_n+\frac{h}{2}k_2)k_4=f(x_n+h,y_n+hk_3)y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)其中，k_i为中间变量。
三、常微分方程组的数值解法1.欧拉法对于一个二阶常微分方程组：\begin{cases} y'_1=f_1(x,y_1,y_2) \\ y'_2=f_2(x,y_1,y_2)\end{cases}其欧拉公式为：\begin{cases} y_{n+1,1}=y_{n,1}+hf_1(x_n,y_{n,1},y_{n,2}) \\y_{n+1,2}=y_{n,2}+hf_2(x_n,y_{n,1},y_{n,2}) \end{cases}其中，x_n和y_{n,i}(i=1, 2)为第n个时间点上x和y_i的取值。
常微分方程的数值解法
常微分方程的数值解法常微分方程是研究变量的变化率与其当前状态之间的关系的数学分支。
它在物理、工程、经济等领域有着广泛的应用。
解常微分方程的精确解往往十分困难甚至不可得，因此数值解法在实际问题中起到了重要的作用。
本文将介绍常见的常微分方程的数值解法，并比较其优缺点。
1. 欧拉方法欧拉方法是最简单的数值解法之一。
它基于近似替代的思想，将微分方程中的导数用差商近似表示。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）根据微分方程的定义使用近似来计算下一个点的值。
欧拉方法的计算简单，但是由于误差累积，精度较低。
2. 改进欧拉方法为了提高欧拉方法的精度，改进欧拉方法应运而生。
改进欧拉方法通过使用两个点的斜率的平均值来计算下一个点的值。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）根据微分方程的定义使用近似来计算下一个点的值。
改进欧拉方法相较于欧拉方法而言，精度更高。
3. 龙格-库塔法龙格-库塔法（Runge-Kutta）是常微分方程数值解法中最常用的方法之一。
它通过迭代逼近精确解，并在每一步中计算出多个斜率的加权平均值。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）计算各阶导数的导数值。
（4）根据权重系数计算下一个点的值。
与欧拉方法和改进欧拉方法相比，龙格-库塔法的精度更高，但计算量也更大。
4. 亚当斯法亚当斯法（Adams）是一种多步法，它利用之前的解来近似下一个点的值。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）通过隐式或显式的方式计算下一个点的值。
亚当斯法可以提高精度，并且比龙格-库塔法更加高效。
5. 多步法和多级法除了亚当斯法，还有其他的多步法和多级法可以用于解常微分方程。
多步法通过利用多个点的值来逼近解，从而提高精度。
而多级法则将步长进行分割，分别计算每个子问题的解，再进行组合得到整体解。
常微分方程的数值解法全文
第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 <<.。
常微分方程数值解法
第八章 常微分方程的数值解法一．内容要点考虑一阶常微分方程初值问题:⎪⎩⎪⎨⎧==00)(),(y x y y x f dx dy微分方程的数值解：设微分方程的解y (x )的存在区间是[a,b ]，在[a,b ]内取一系列节点a= x 0< x 1<…< x n =b ,其中h k =x k+1-x k ；(一般采用等距节点，h=(b-a)/n 称为步长)。
在每个节点x k 求解函数y(x)的近似值:y k ≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。
用数值方法，求得f(x k )的近似值y k ，再用插值或拟合方法就求得y(x)的近似函数。
(一)常微分方程处置问题解得存在唯一性定理对于常微分方程初值问题:⎪⎩⎪⎨⎧==00)(),(y x y y x f dx dy如果：(1) 在B y y A x x 00≤-≤≤,的矩形内),(y x f 是一个二元连续函数。
(2) ),(y x f 对于y 满足利普希茨条件，即2121y y L y x f y x f -≤-),(),(则在C x x 0≤≤上方程⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy的解存在且唯一，这里C=min((A-x 0),x 0+B/L),L 是利普希茨常数。
定义：任何一个一步方法可以写为),,(h y x h y y k k k 1k Φ+=+，其中),,(h y x k k Φ称为算法的增量函数。
收敛性定理：若一步方法满足： (1)是p 解的.(2) 增量函数),,(h y x k k Φ对于y 满足利普希茨条件.(3) 初始值y 0是精确的。
则),()()(p h O x y kh y =-kh =x -x 0,也就是有0x y y lim k x x kh 0h 0=--=→)((一)、主要算法 1．局部截断误差局部截断误差：当y(x k )是精确解时，由y(x k )按照数值方法计算出来的1~+k y 的误差y (x k+1)- 1~+k y 称为局部截断误差。
计算方法 第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
用欧拉方法求解常微分方程组
用欧拉方法求解常微分方程组欧拉方法是求解常微分方程组的一种数值解法。
当解析解难以求得时，我们可以通过数值方法来近似求解。
欧拉方法是一种简单粗暴的方法，它的基本思想是将微分方程在一段区间内用直线逼近，然后通过不断迭代来逼近解。
我们来看一个简单的一阶常微分方程，y’=y，初始条件为y(0)=1。
通过欧拉方法，我们可以将其离散化为y(n+1)=y(n)+hf(y(n))，其中h为步长。
我们从初始点开始，以步长h不断向前迭代，即可得到y(1)=y(0)+hy(0)=2，y(2)=y(1)+hy(1)=2.718，y(3)=y(2)+hy(2)=4.669…，一次迭代可以得到一个新的解点，通过不断迭代，我们可以逼近真正的解。
欧拉方法的原理简单易懂，但是由于直线逼近的误差会随着步长的增大而增大，所以我们需要选择合适的步长来保证精度。
此外，欧拉方法只适用于简单的一阶线性常微分方程，对于高阶、非线性方程的求解并不是很有效。
接下来，我们来看如何用欧拉方法求解二阶常微分方程。
我们将其离散化，得到y(n+1)=y(n)+hy’(n)，y’(n+1)=y’(n)+hf(y’(n),y(n))，其中y(n)表示y在第n个步长的值，y’(n)表示y’在第n个步长的值，f(y’,y)表示微分方程右侧的函数，它可以由高阶常微分方程离散化得到。
同样的，通过不断迭代，我们可以得到y(n)和y’(n)的数值解。
不难发现，欧拉方法的复杂度为O(N)，其中N表示步数，所以当步数较大时，欧拉方法的计算复杂度较高，同时误差也会增大。
因此，我们需要选择合适的步长来保证精度，一般可以使用控制误差的算法来自适应地调整步长。
总的来说，欧拉方法是一种简单易懂，实现容易的数值解法。
虽然它的精度在高精度和高维度的情况下并不是很高，但对于一些简单的常微分方程组求解还是非常有效的。
常微分方程数值解法
常微分方程数值解法常微分方程是研究函数的导数与自变量之间的关系的数学分支，广泛应用于物理、工程、生物等领域的建模与分析。
在实际问题中，我们常常遇到无法通过解析方法求得精确解的常微分方程，因此需要利用数值解法进行求解。
本文将介绍几种常用的常微分方程数值解法。
一、欧拉方法(Euler's Method)欧拉方法是最基本的数值解法之一。
它的思想是将微分方程转化为差分方程，通过逐步逼近解的方式求得数值解。
具体步骤如下：1. 将微分方程转化为差分方程：根据微分方程的定义，可以得到差分方程形式。
2. 选择步长：将自变量范围进行离散化，确定步长h。
3. 迭代计算：根据差分方程递推公式，利用前一步的数值解计算后一步的数值解。
二、改进的欧拉方法(Improved Euler's Method)改进的欧拉方法通过使用欧拉方法中的斜率来进行更准确的数值计算。
具体步骤如下：1. 计算欧拉方法的斜率：根据当前节点的数值解计算斜率。
2. 根据斜率计算改进的数值解：将得到的斜率代入欧拉方法的递推公式中，计算改进的数值解。
三、龙格-库塔方法(Runge-Kutta Method)龙格-库塔方法是一类常微分方程数值解法，其中最著名的是四阶龙格-库塔方法。
它通过计算各阶导数的加权平均值来逼近解，在精度和稳定性方面相对较高。
具体步骤如下：1. 计算每一步的斜率：根据当前节点的数值解计算每一步的斜率。
2. 计算权重：根据斜率计算各个权重。
3. 计算下一步的数值解：根据计算得到的权重，将其代入龙格-库塔方法的递推公式中，计算下一步的数值解。
四、多步法(多步差分法)多步法是需要利用多个前面节点的数值解来计算当前节点的数值解的数值方法。
常见的多步法有Adams-Bashforth法和Adams-Moulton法。
具体步骤如下：1. 选择初始值：根据差分方程的初始条件，确定初始值。
2. 迭代计算：根据递推公式，利用前面节点的数值解计算当前节点的数值解。
- 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];