当前位置：文档之家 › MATlAB 数值求解ODE方程

MATlAB 数值求解ODE方程

  • 格式：pdf
  • 大小：292.07 KB
  • 文档页数：10

下载文档原格式

  / 10

合集下载

重要：MATLAB常微分方程(组)数值解法

重要：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];

matlab解二阶偏微分方程

matlab解二阶偏微分方程

<div align="center">一、Matlab求解二阶偏微分方程(ODE)的基本步骤</div>1. 数学模型：首先要确定求解的方程是哪一类的偏微分方程(ODE)，然后建立其对应的数学模型，使其符合这类微分方程的形式;2. 确定边界条件：确定迭代范围$[a,b]$,边界条件函数 $y(a)=\alpha$ 、$y(b)=\beta$;3. 写出Matlab程序：在该类ODE中，通常会有某一种常用的数值求解方法，一般使用微分方程求解器(ODE)，如ode45等;4. 获得实际结果：开始编写Matlab程序，完成参差和参数的输入以后，可以运行Matlab程序，然后求得结果，再用图像表示出来。

<div align="center">二、具体求解</div>$$\frac{d^2y}{dx^2}+y=6sin(2x)$$微分方程为二阶常微分方程，求解条件如下：$[a,b]=[0,\pi], y(0)=1,y(\pi)=3.$（1）Matlab函数表达式首先建立与二阶非齐次线性常微分方程相符合的数学模型，其Matlab函数表达式为$$ f(x,y,y')=\frac{dy}{dx}-y'-6sin2x $$其中，$y=y(x)$；（2）函数程序在Matlab中，定义函数程序 $myode.m$ ，此程序返回右端函数 $f(x,y,y')$ 的值表达式，程序内容如下。

```MATLAB% 右端函数程序function dy=myode(x,y)dy=[y(2);-y(2)-6*sin(2*x)];end```（3）调用Matlab函数olvede45调用Matlab函数 solvede45 求解二阶ODN，程序内容如下：```MATLAB% 主程序求解% maxstep表示分裂的步长大小% Tolerence表示误差，控制求解精度Maxstep=0.25;Tolerence=1e-4;a=0;b=pi;y0=[1;0];[x,y] = ode45('myode',[a,b],y0,options);```（4）结果展示输入参数之后，运行Matlab程序，得到如下图：![](../images/matlab_2_diff.png)此图为$y(x)$随$x$变化的曲线，可以看出，二阶偏微分的求解结果满足了边界条件，即$y(0)=1,y(\pi)=3$ ，如图中红色圆点所示。

ODE问题解析解及数值解的matlab实现

ODE问题解析解及数值解的matlab实现

龙源期刊网 ODE问题解析解及数值解的matlab实现作者：于丽妮来源：《电脑知识与技术》2012年第14期利用数学手段研究自然现象和社会现象或解决工程技术问题，如揭示质点运动规律的Newton第二定律和刻画回路电流或电压变化规律的基尔霍夫回路定律等。

一般先要建立数学模型，再对数学模型进行简化和求解，然后结合实际问题对结果进行分析和讨论。

数学模型中最常见的表达方式，是包含自变量和未知函数以及未知函数导数的的函数方程。

即微分方程，通常包括常微分方程（ODE）和偏微分方程（PDE）。

ODE是常微分方程的英文缩写，即ordinary diffrential equation，如果在微分方程中，自变量的个数只有一个，这就是ODE方程，例如形如F(x,y,y’,y")=0的方程就是一个二阶ODE方程。

对于常微分方程来说，能够求出解析解的是有限的，即不能通过等式左右同时积分求不定积分的方法求出，特别是高阶方程和偏微分方程，有的解还需要用插值方法继续计算，因此用解析法来求微分方程的数值解有时是不适宜的。

但是在给定条件下，可以求出常微分方程的数值解来逼近数值解，以完成对模型的研究。

这就要求我们必须研究微分方程的解法，既要研究微分方程的解析解法（精确解），更要研究微分方程的数值解法（近似解）。

ODE问题通常是对一动态过程（动态系统、动力系统）演化规律的描述，求解ODE问题就是要了解和掌握动态过程演化规律。

本文主要研究ODE问题解析解及数值解的matlab实现。

微分方程是一门实用性很强的学科，对于理论研究以及生活实际提出的微分方程问题，能够求出解析解的方程是有限的，因此将计算机应用到求微分方程解的问题是必要的。

本文详细介绍了常微分方程中解析解及数值解解法的理论基础，并用matlab加以实现。

matlab中ode23的使用示例

matlab中ode23的使用示例

在MATLAB中，ode23函数是用于求解常微分方程（ODE）的数值解法之一。

它采用二阶或三阶Runge-Kutta算法进行求解。

以下是ode23函数的使用示例：
假设我们要求解以下常微分方程：
dy/dt = y - t^2 + 1
初始条件是y(0) = 0.5，求解区间为[0, 2]。

在MATLAB中，首先需要编写一个函数来描述这个常微分方程。

可以创建一个名为"odeFunc.m"的函数文件，包含以下内容：
matlab
function dy = odeFunc(t, y)
dy = y - t^2 + 1;
end
然后，在MATLAB命令窗口或脚本文件中，使用ode23函数来求解该常微分方程。

示例代码如下：
matlab
% 定义求解区间、初始条件和ODE函数
tspan = [0, 2];
y0 = 0.5;
odeFunc = @odeFunc;
% 使用ode23函数求解ODE
[t, y] = ode23(odeFunc, tspan, y0);
% 绘制结果图形
plot(t, y);
xlabel('Time (t)');
ylabel('Solution (y)');
title('Solution of ODE using ode23');
运行上述代码后，MATLAB将会计算并返回在指定求解区间内的常微分方程数值解，并将结果绘制成图形。

图形中的x轴表示时间（t），y轴表示解（y）。

matlab_常微分方程数值解法

matlab_常微分方程数值解法
d2x 2x2 0
dt 2
简朴问题可以求得解析解，多数实际问题靠数值求解 。
第4页
一阶常微分方程(ODE )初值问题 ： ODE ：Ordinary Differential Equation
dy
f
(x,
y)
dx
x0 x xn
y(x0 ) y0
数值解法就是求y(x)在某些分立旳节点 xn 上旳近似值 yn，用以近似y(xn)
x0
y0
x1 f y(x), x dx
x0
x2 f y(x), x dx
x1
y(x1) f y(x1), x1 h
第17页
同样，在[x0,xn+1] ，积分采用矩形近似，得：
y(xn1) y0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
yn y(xn )
第5页
2、欧拉近似办法
2.1 简朴欧拉(L.Euler, 1707-1783)办法。
dy
dx
f
(y, x)
y(x0 ) y0
欧拉数值算法就是由初值通过递推求解，递推求解
就是从初值开始，后一种函数值由前一种函数值得到。核 心是构造递推公式。
y0 y1 y2 yn
第6页
i 1,2,...
第36页
没有一种算法可以有效地解决所有旳 ODE 问题，因此 MATLAB 提供了多种ODE函数。
函数 ODE类
特点
阐明
型
ode45
非刚性 单步法；4，5 阶 R-K 措施；合计 大部分场合旳首选措施
截断误差为 (△x)3
ode23
非刚性 单步法；2，3 阶 R-K 措施；合计 使用于精度较低旳情形

matlab数值求解常微分方程快速方法

matlab数值求解常微分方程快速方法

MATLAB是一种用于科学计算和工程应用的高级编程语言和交互式环境。

它在数学建模、模拟和分析等方面有着广泛的应用。

在MATLAB 中，常微分方程的数值求解是一个常见的应用场景。

在实际工程问题中，通常需要对常微分方程进行数值求解来模拟系统的动态行为。

本文将介绍MATLAB中对常微分方程进行数值求解的快速方法。

1. 基本概念在MATLAB中，可以使用ode45函数来对常微分方程进行数值求解。

ode45是一种常用的Runge-Kutta法，它可以自适应地选取步长，并且具有较高的数值精度。

使用ode45函数可以方便地对各种类型的常微分方程进行求解，包括一阶、高阶、常系数和变系数的微分方程。

2. 函数调用要使用ode45函数进行常微分方程的数值求解，需要按照以下格式进行函数调用：[t, y] = ode45(odefun, tspan, y0)其中，odefun表示用于描述微分方程的函数，tspan表示求解的时间跨度，y0表示初值条件，t和y分别表示求解得到的时间序列和对应的解向量。

3. 示例演示为了更好地理解如何使用ode45函数进行常微分方程的数值求解，下面我们以一个具体的例子来进行演示。

考虑如下的一阶常微分方程：dy/dt = -2*y其中，y(0) = 1。

我们可以编写一个描述微分方程的函数odefun：function dydt = odefun(t, y)dydt = -2*y;按照上述的函数调用格式，使用ode45函数进行求解：tspan = [0 10];y0 = 1;[t, y] = ode45(odefun, tspan, y0);绘制出解曲线：plot(t, y);4. 高级用法除了基本的函数调用方式外，MATLAB中还提供了更多高级的方法来对常微分方程进行数值求解。

可以通过设定选项参数来控制数值求解的精度和稳定性，并且还可以对刚性微分方程进行求解。

5. 性能优化在实际工程应用中，常常需要对大规模的常微分方程进行数值求解。

matlab求微分方程数值解

matlab求微分方程数值解利用matlab求微分方程数值解是一种常用的数学计算方法。

在实际工程和科学研究中，许多问题都可以用微分方程来描述，但是解析解往往难以求得，因此需要用数值方法求解微分方程。

求解微分方程的数值方法有很多种，其中比较常用的是欧拉法和龙格-库塔法。

欧拉法是一种基本的数值方法，它采用离散化的方法将微分方程转化为差分方程，然后通过迭代来求出数值解。

欧拉法的具体步骤是：首先将自变量和因变量离散化，然后利用微分方程的定义式将微分方程转化为差分方程，最后通过迭代求出数值解。

欧拉法的优点是简单易懂，但是精度较低，容易产生误差。

龙格-库塔法是一种高阶数值方法，它将微分方程转化为一系列的差分方程，并采用递推的方法求解数值解。

龙格-库塔法的优点是精度高，收敛速度快，适用于求解复杂的微分方程。

但是龙格-库塔法的计算量较大，需要进行多次计算，计算时间较长。

在使用matlab求解微分方程时，可以直接调用matlab中的ode 函数来求解微分方程。

ode函数是matlab中内置的求解微分方程的函数，它支持多种数值方法，包括欧拉法和龙格-库塔法等。

使用ode函数可以简化求解微分方程的过程，提高计算效率。

在使用ode函数求解微分方程时，需要先定义微分方程的函数表达式，然后将函数表达式作为参数传入ode函数中。

ode函数会自动选择合适的数值方法来求解微分方程，并返回数值解。

通过调整ode函数的参数，可以进一步提高求解微分方程的精度和计算效率。

除了ode函数外，matlab中还有很多其他的数值计算函数，如dsolve函数、pdepe函数等，它们可以用来求解不同类型的微分方程。

在实际应用中，需要根据具体问题选择合适的数值方法和函数来求解微分方程。

利用matlab求解微分方程数值解是一种常用的数学计算方法，可以通过调用matlab中的内置函数来实现。

在选择数值方法和函数时需要考虑精度和计算效率等因素，以便更好地解决实际问题。

Matlab ode函数 微分方程的数值解

[T,Y] = ode45(odefun,tspan,y0)
odefun是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名
tspan是区间[t0 tf]或者一系列散点[t0,t1,...,tf]
y0是初始值向量《Simulink与信号处理》
T返回列向量的时间点
Y返回对应T的求解列向量
[T,Y] = ode45(odefun,tspan,y0,options)
options是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等
[T,Y,TE,YE,IE] =ode45(odefun,tspan,y0,options)
在设置了事件参数后的对应输出
TE事件发生时间
YE事件解决时间
IE The index i of the event functionthat vanishes.
sol =ode45(odefun,[t0 tf],y0...)
sol结构体输出结果
编辑本段
示例
求解一阶常微分方程
odefun=@(t,y) (y+3*t)/t^2; %定义函数
tspan=[1 4]; %求解区间
y0=-2; %初值
[t,y]=ode45(odefun,tspan,y0);plot(t,y) %作图
title('t^2y''=y+3t,y(1)=-2,1<t<4') legend('t^2y''=y+3t')
xlabel('t')
ylabel('y')=y+3*t','y(1)=-2')

matlab迭龙格库塔法解常微分方程

一、介绍迭龙格-库塔法（Runge-Kutta method）是一种数值求解常微分方程（ODE）的常用方法。

它是由卡尔·迭龙格（Carl Runge）和马丁·威尔黑尔姆·库塔（Wilhelm Kutta）在20世纪初提出的，该方法以两位数值分析家的名字来命名。

二、简单描述迭龙格-库塔法是通过数值逼近的方式，来计算常微分方程的近似解。

它是一种显式求解方法，适用于解非线性常微分方程和具有较大阶数的常微分方程。

三、数学原理迭龙格-库塔法主要是通过将微分方程转化为差分方程，利用数值解的方式来逼近微分方程的解。

它是一种显式方法，通过不断迭代得到下一个时间步的近似解。

四、matlab中的应用在matlab中，可以使用ode45函数来调用迭龙格-库塔法求解常微分方程。

ode45函数是matlab中集成的一个函数，通过调用ode45函数，可以直接求解常微分方程的数值解。

五、实例演示下面通过一个简单的例子来演示如何使用matlab中的ode45函数来求解常微分方程。

我们考虑一个简单的一阶常微分方程：dy/dt = -y初始条件为y(0) = 1。

在matlab中，可以通过以下代码来求解该微分方程：```定义微分方程的函数function dydt = myode(t, y)dydt = -y;调用ode45函数求解[t, y] = ode45(myode, [0, 5], 1);plot(t, y);```运行以上代码，即可得到微分方程的数值解，并通过绘图来展示解的变化。

六、总结迭龙格-库塔法是一种常用的数值解常微分方程的方法，它在matlab中有较为方便的调用方式。

通过ode45函数，可以快速求解常微分方程的数值解，并通过绘图来展示结果。

希望本篇文章对读者有所帮助，谢谢阅读。

七、应用场景和优势在实际应用中，迭龙格-库塔法广泛应用于各种科学和工程领域，如物理学、化学、生物学、经济学等。

matlab第6章ode


1 L 0
u
i
状态方程以两个状态元件 i和uo作状态变量 21
+S
i
ui
D
-
dT
L
+
C u0 R
-
0
L
di dt
u0
i
C
d u0
dt
u0
R
T
(2)当S=off：时间长度为: （1-d)*T:
di dt
u0
L
du0
dt
i C
u0
RC
22
S=OFF
u
•
i
•
0
0
1 C
1
1. ode23
在MATLAB中，函数ode23采用2-3阶龙格-库塔法 求解微分方程。
[t,y]=ode23(odefun,tspan,y0) [t,y]=ode23(odefun,tspan,y0,options) odefun：定义微分方程的形式y’=f(t,y) tspan=[t0,tfinal]：表示微分方程的积分限从t0（始值） 到tfinal（终值），该积分限也可以是一些离散的点。 y0:初始状态列向量 options：积分参数，包括‘RelTol’（相对误差）和 ‘AbsTol’（绝对误差），可省略。
IR Vi
L C UC
12
（1）分析：根据电路分析，可以得出微分方程
RI(t) VL(t) VC (t) Vi (t)
I(t) C dVC (t) dt
VL (t)
L
dI(t) dt
LC
d
2VC (t) dt 2
dI(t) L dt RI(t) VC (t) Vi (t)
LC
  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

pk = xk + hf (tk, xk), tk+1 = tk + h, h
xk+1 = xk + 2 [f (tk, xk) + f (tk+1, pk)] .
4
1.3 Matlab Builtin ODE Solvers
In addition there are many other methods for approximating solutions to ordinary differential equations, but due to a lack of time left in the semester I will just introduce you to Matlabs built-in Runge-Kutta solver ode45 and show you how it works.
s1 = f (tk, xk)
h
h
s2 = f tk + 2 , xk + 2 s1
xk+1 = xk + hs2
tk&s obtained as follows: First we replace the approximation x(t + h) − x(t)
next section. They are all of higher order than Euler’s method.
1.2 Modified Euler’s Methods
Once again we consider (1) and apply the fundamental theorem of calculus to get
t2
t2
f (t, x(t)) dt = x (t) dt = x(t2) − x(t1).
t1
t1
This implies that
t2
x(t2) = x(t1) + f (t, x(t)) dt.
t1
Now we consider two possible ways to approximate the integral on the right above.
x (t) = e−1/2 t2 .
We build the m-file eul1d.m and save it to disk (use cd T:\ to change to the T drive)
% Euler’s method 1 clear all defn=inline(’-t*x’,’t’,’x’); a=0; b=5; N=100; xa=1; h=(b-a)/N; t=a+h*(0:N); LT=length(t); x(1)=xa;
−tx(t)
x (t) =
, x(0) = 1, x (0) = 1, on 0 < t < 5.
(x (t)2 + 1)
Now we build an m-file eul2d.m to solve the problem.
clear % Euler’s method 2d clear defn=inline(’[y;-t*x/(y^2+1)]’,’t’,’x’,’y’); a=0; b=5; N=200; h=(b-a)/N; t=a+h*(0:N); LT=length(t); xa=1; xpa=1; w(:,1)=[xa;xpa];
x = f (t, x), t ∈ (a, b)
(1)
x(a) = xa
(2)
we observe that defining tj = a + hj, for j = 0, 1, 2, · · · , N where h = (b − a)/N we have
x(tj+1) − h
x(tj )
≈
f (tj, x(tj)).
for j=2:LT w(:,j)=w(:,j-1)+h*defn(t(j-1),w(1,j-1),w(2,j-1));
end x=w(1,:); plot(t,x,’LineWidth’,2)
grid on 2
We plot the result of executing this Matlab code in the following figure.
x = f (t, x, x ), t ∈ [a, b]
x(a) = xa x (a) = xpa
If we set y = x , then the sxstem can be written as
d dt
x y
=
y f (t, x, y)
x(a) = xa
y(a) = xpa
The numerical approximation can now be carrried out just as above. As we have already mentioned it is very difficult to find exact solutions to most interesting problems. The next example is a second order initial value problem. Consider
Math 4330 Sec. 1, Matlab Assignment # 4 , April 26, 2006 Name
1 Numerical Solution of ODEs Using Matlab
1.1 Euler’s Method
Euler’s one step method is undoubtedly the simplest method for approximating the solution to an ordinary differential equation. It goes something like this: Given a first order initial value problem
The basic syntax for using ode45 is the following
[t,x] = ode45(’xdot’,tspan,x0)
where xdot.m is a function file containing the the right hand side of the differential equation, tspan is a vector of time values at which you want to obtain the solution and x0 is the initial condition.
x (t) = h
3
by the more accurate approximation (see the figure)
h x(t + h) − x(t)
x t+ =
2
h
Then we apply Taylor’s formula of order one
g(ξ) = g(a) + g (a)(ξ − a) + O(((ξ − a)2),
xj+1 = xj + f (tj, xj), j = 0, 2, · · · , (N − 1).
As an example let us consider the IVP
x (t) = −tx(t), x(0) = 1 on 0 < t < 5.
The exact solution is given by
The Midpoint Rule:
The midpoint method uses Euler to step halfway across the interval, evaluates the function at this
intermediate point, then uses that slope to take the actual step.
2
2
2
h
h
h
x(t + h) ≈ x(t) + f t + , x(t) + f (t, x(t)) .
2
2
2
The Trapezoid Rule: For the trapezoid rule we approximate the integral with step size h = t2 − t1 using
h x(t2) ≈ x(t1) + 2 [f (t1, x(t1)) + f (t2, x(t2))].
Unfortunately, we do not know x(t2) on the right side, so we use, for example, Euler’s method to approximate it: x(t2) ≈ x(t1) + hf (t1, x(t1)). We thus obtain a numerical procedure by denoting x(tk) by yk for each k and setting
for j=2:LT x(j)=x(j-1)+h*defn(t(j-1),x(j-1));
end xact_sol=exp(-t.^2/2); plot(t,x,’b’,t,xact_sol,’r--’,’LineWidth’,2)
Execute the file in the Workspace by typing the first name of the file without the .m extension. We plot the result of executing this Matlab code, on 0 ≤ t ≤ 5 with initial condition x(0) = 1, in the following figure.