1. 全栈程序员必看首页

matlab中ode45函数解二阶微分方程_matlab求常微分方程组

全栈程序员-站长 • • 未分类 • 阅读 4

matlab中ode45函数解二阶微分方程_matlab求常微分方程组Matlab微分方程求解并绘制曲线1.用dsolve()求解>>clear>>clc>>symsy(t)>>Dy=diff(y,1)Dy(t)=diff(y(t),t)>>y=dsolve(Dy==y-2*t/y,y(0)==1)y=(2*t+1)^(1/2)>>t=0:0.1:4;>>y=eval(y);>>plot(t,

大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全家桶1年46，售后保障稳定

Matlab 微分方程 ode45 求解并绘制曲线

2. 用 ode45() 求解

2.1 ode45() 函数用法

[t, Xt] = ode45(odefun, tspan, X0)

  • odefun 是函数句柄，可以是函数文件名，匿名函数句柄或内联函数名
  • tspan 是区间 [t0 tfinal] 或者一系列散点[t0,t1,…,tf]
  • X0 是初始值向量
  • t 返回列向量的时间点
  • Xt 返回对应T的求解列向量

2.2 示例：求解一阶微分方程

求解单变量微分方程的解
x ˙ ( t ) = 2 ∗ x ( t ) \dot{x}(t) = 2 * x(t) x˙(t)=2∗x(t)

2.2.1 Matlab 代码如下

clear;
clc;

% 程序主函数代码如下：
t0 = 0;
tfinal = 10;
X0 = 1;

[t, Xt] = ode45(@SunFun, [t0 tfinal], X0);

% 绘制结果图
plot(t,Xt)
grid

% 微分方程函数，状态导数
function xdot = SunFun(t,x)
% 导数关系式
xdot = 2 * x;

end

Jetbrains全家桶1年46，售后保障稳定

2.2.2 代码效果

在这里插入图片描述

2.3 示例：求解矩阵一阶微分方程

求解一阶微分方程
X ˙ ( t ) = − L ∗ X ( t ) \dot{X}(t) = -L * X(t) X˙(t)=−L∗X(t)

其中， X ( t ) X(t) X(t) 是列向量。

2.3.1 Matlab 代码

clear;
clc;

% 程序主函数代码如下：
t0 = 0;
tfinal = 10;
X0 = [0; 1; 3; 3; 5; 6;];

[t, Xt] = ode45(@SunFun, [t0 tfinal], X0);

% 绘制结果图
plot(t,Xt(:,1), t,Xt(:,2), t,Xt(:,3), t,Xt(:,4), t,Xt(:,5), t,Xt(:,6))
grid

% 微分方程函数，状态导数
function xdot = SunFun(t,x)

D = [1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1;];
A = [0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1; 1 0 0 0 0 0;];
L = D - A;

% 导数关系式
xdot = -L * x;

end

2.3.2 代码效果

在这里插入图片描述

2.4 示例：将上述示例代码写成两个函数

2.4.1 主函数如下

clear;
clc;

% 程序主函数代码如下：
t0 = 0;
tfinal = 10;
X0 = [0; 1; 3; 3; 5; 6;];

[t, Xt] = ode45('SunFun', [t0 tfinal], X0);

% 绘制结果图
plot(t,Xt)
grid

2.4.2 子函数如下

% 微分方程函数，状态导数
function xdot = SunFun(t,x)

D = [1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1;];
A = [0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1; 1 0 0 0 0 0;];
L = D - A;

% 导数关系式
xdot = -L * x;

From: P4 控制系统数学模型-《Matlab/Simulink与控制系统仿真》程序指令总结

Ref: 【MATLAB】关于ode45的一部分用法的函数编写方式

1. ode45-官方释义

1.1 语法 / 说明

[t,y] = ode45(odefun,tspan,y0)

[t,y] = ode45(odefun,tspan,y0)（其中 tspan = [t0 tf]）求微分方程组 y′=f(t,y) 从 t0 到 tf 的积分，初始条件为 y0。解数组 y 中的每一行都与列向量 t 中返回的值相对应。

所有 MATLAB® ODE 求解器都可以解算 y′=f(t,y) 形式的方程组，或涉及质量矩阵 M(t,y)y′=f(t,y) 的问题。求解器都使用类似的语法。ode23s 求解器只能解算质量矩阵为常量的问题。ode15s 和 ode23t 可以解算具有奇异质量矩阵的问题，称为微分代数方程 (DAE)。使用 odeset 的 Mass 选项指定质量矩阵。

ode45 是一个通用型 ODE 求解器，是您解算大多数问题时的首选。但是，对于刚性问题或需要较高准确性的问题，其他 ODE 求解器可能更适合。有关详细信息，请参阅选择 ODE 求解器。

[t,y] = ode45(odefun,tspan,y0,options)

[t,y] = ode45(odefun,tspan,y0,options) 还使用由 options（使用 odeset 函数创建的参数）定义的积分设置。例如，使用 AbsTol 和 RelTol 选项指定绝对误差容限和相对误差容限，或者使用 Mass 选项提供质量矩阵。

[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)
[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options) 还求 (t,y) 的函数（称为事件函数）在何处为零。在输出中，te 是事件的时间，ye 是事件发生时的解，ie 是触发的事件的索引。

对于每个事件函数，应指定积分是否在零点处终止以及过零方向是否重要。为此，请将 ‘Events’ 属性设置为函数（例如 myEventFcn 或 @myEventFcn），并创建一个对应的函数：[value,isterminal,direction] = myEventFcn(t,y)。有关详细信息，请参阅 ODE 事件位置。

sol = ode45(___)

sol = ode45(___) 返回一个结构体，您可以将该结构体与 deval 结合使用来计算区间 [t0 tf] 中任意点位置的解。您可以使用上述语法中的任何输入参数组合。

1.2 示例

1.2.1 具有一个解分量的 ODE

在对求解器的调用中，可将只有一个解分量的简单 ODE 指定为匿名函数。该匿名函数必须同时接受两个输入 (t,y)，即使其中一个输入未使用也是如此。

解算 ODE
y ′ = 2 t y’ = 2t y′=2t

使用时间区间 [0,5] 和初始条件 y0 = 0。

tspan = [0 5];
y0 = 0;
[t,y] = ode45(@(t,y) 2*t, tspan, y0);


% 对解绘图。
plot(t,y,'-o')

在这里插入图片描述

1.2.2 van der Pol 方程为二阶 ODE

y 1 ′ ′ − μ ( 1 − y 1 2 ) y 1 ′ + y 1 = 0 y”_1 – \mu (1-y_1^2)y’_1 + y_1 = 0 y1′′​−μ(1−y12​)y1′​+y1​=0

其中， μ > 0 \mu>0 μ>0 为标量参数。通过执行 y 1 ′ = y 2 y’_1=y_2 y1′​=y2​ 代换，将此方程重写为一阶 ODE 方程组。生成的一阶 ODE 方程组为
y 1 ′ = y 2 y 2 ′ = μ ( 1 − y 1 2 ) y 2 − y 1 y’_1 = y_2 \\ y’_2 = \mu(1-y_1^2)y_2-y_1 y1′​=y2​y2′​=μ(1−y12​)y2​−y1​

函数文件 vdp1.m 代表使用 μ = 1 \mu = 1 μ=1 的 van der Pol 方程。变量 y 1 y_1 y1​ 和 y 2 y_2 y2​ 是二元素向量 dydt 的项 y ( 1 ) y(1) y(1) 和 y ( 2 ) y(2) y(2)。

function dydt = vdp1(t,y)
	dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
end

使用 ode45 函数、时间区间 [0 20] 和初始值 [2 0] 来解算该 ODE。生成的输出即为时间点 t t t 的列向量和解数组 y y y。 y y y 中的每一行都与 t t t 的相应行中返回的时间相对应。 y y y 的第一列与 y 1 y_1 y1​ 相对应，第二列与 y 2 y_2 y2​ 相对应。

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

% 绘制 $y_1$ 和 $y_2$ 的解对 t 的图。
plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

整体写成一个函数

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

% 绘制 $y_1$ 和 $y_2$ 的解对 t 的图。
plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

function dydt = vdp1(t,y)
    p = y(1);
    dp = y(2);
	ddp = (1-p^2)*dp-p;
    dydt = [dp; ddp];
end

再做一下转换，换成常见的二阶智能体形式：
p ˙ = v v ˙ = ( 1 − p 2 ) ∗ v − p \dot{p}=v \\ \dot{v}=(1-p^2)*v-p p˙​=vv˙=(1−p2)∗v−p

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

% 绘制 $y_1$ 和 $y_2$ 的解对 t 的图。
plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

function out = vdp1(t,y)
    p = y(1,:);
    v = y(2,:);
	ddp = (1-p^2)*v-p;
    out = [v; ddp];
end

在这里插入图片描述

因为对ode的使用方法（求二阶微分方程）还是不够熟悉，因此，再做一次离散化处理，来验证下自己的离散化方法是否正确。

在这里插入图片描述

图像变化趋势没有错误，证明离散化的思路正确，再把时间间隔由 d T = 0.1 dT = 0.1 dT=0.1 调整为 d T = 0.01 dT=0.01 dT=0.01，结果如下
在这里插入图片描述

最后，附上完整版代码：

tBegin = 0;
tFinal = 20;
dT = 0.01;
times = (tFinal - tBegin)/dT;

p(:,1) = 2;
v(:,1) = 0;
t(:,1) = 0;

for i = 1:times
    t(i+1,1) = t(i,1) + dT;
    p(i+1,1) = p(i,1) + dT * v(i,1);
    v(i+1,1) = v(i,1) + dT * ((1-p(i,1)^2)*v(i,1)-p(i,1));
end

y = [p,v];

% 绘制 $y_1$ 和 $y_2$ 的解对 t 的图。
plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

1.2.3 向 ODE 函数传递额外的参数

ode45 仅适用于使用两个输入参数（ t t t 和 y y y）的函数。但是，通过在函数外部定义参数并在指定函数句柄时传递这些参数，可以传入额外参数。

解算 ODE
y ′ ′ = A B t y y” = \frac{A}{B}ty y′′=BA​ty

将该方程重写为一阶方程组可以得到
y 1 ′ = y 2 y 2 ′ = A B   t   y 1 y’_1 = y_2 \\ y’_2 = \frac{A}{B}\ t\ y_1 y1′​=y2​y2′​=BA​ t y1​

odefcn.m 将此方程组表示为接受四个输入参数（t、y、A 和 B）的函数。

function dydt = odefcn(t,y,A,B)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = (A/B)*t.*y(1);

使用 ode45 解算 ODE。指定函数句柄，使其将 A 和 B 的预定义值传递给 odefcn。

A = 1;
B = 2;
tspan = [0 5];
y0 = [0 0.01];
[t,y] = ode45(@(t,y) odefcn(t,y,A,B), tspan, y0);


% 绘制结果。
plot(t,y(:,1),'-o',t,y(:,2),'-.')

在这里插入图片描述

1.3.4 带有时变项的 ODE

请考虑以下带有时变参数的 ODE：

y ′ ( t ) + f ( t ) y ( t ) = g ( t ) . y'(t) + f(t)y(t) = g(t). y′(t)+f(t)y(t)=g(t).

初始条件为 y 0 = 1 y_0 = 1 y0​=1。函数 f(t) 由在时间 ft 时计算的 n×1 向量 f 定义。函数 g(t) 由在时间 gt 时计算的 m×1 向量 g 定义。

创建向量 f 和 g。

ft = linspace(0,5,25);
f = ft.^2 - ft - 3;

gt = linspace(1,6,25);
g = 3*sin(gt-0.25);

编写名为 myode 的函数，该函数通过对 f 和 g 进行插值获取时变项在指定时间的值。将函数保存到您当前的文件夹中，以运行示例的其余部分。

myode 函数接受额外的输入参数以计算每个时间步的 ODE，但 ode45 只使用前两个输入参数 t 和 y。

function dydt = myode(t,y,ft,f,gt,g)
f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t
g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t
dydt = -f.*y + g; % Evaluate ODE at time t

使用 ode45 计算方程在时间区间 [1 5] 内的解。使用函数句柄指定函数，从而使 ode45 只使用 myode 的前两个输入参数。此外，使用 odeset 放宽误差阈值。

tspan = [1 5];
ic = 1;
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[t,y] = ode45(@(t,y) myode(t,y,ft,f,gt,g), tspan, ic, opts);


% 绘制解 y 对时间点 t 的函数图。
plot(t,y)

在这里插入图片描述

1.2.5 计算和扩展结构体

van der Pol 方程为二阶 ODE
y 1 ′ ′ − μ ( 1 − y 1 2 ) y 1 ′ + y 1 = 0 y”_1 -\mu(1-y_1^2) y’_1 + y_1 = 0 y1′′​−μ(1−y12​)y1′​+y1​=0

使用 ode45 以及 μ=1 解算 van der Pol 方程。函数 vdp1.m 随 MATLAB® 一起提供，用于对方程进行编码。指定单个输出以返回包含解信息（如求解器和计算点）的结构体。

tspan = [0 20];
y0 = [2 0];
sol = ode45(@vdp1,tspan,y0)
sol = struct with fields:
     solver: 'ode45'
    extdata: [1x1 struct]
          x: [1x60 double]
          y: [2x60 double]
      stats: [1x1 struct]
      idata: [1x1 struct]

使用 linspace 在区间 [0 20] 内生成 250 个点。使用 deval 计算在这些点上的解。

x = linspace(0,20,250);
y = deval(sol,x);


% 绘制解的第一个分量。
plot(x,y(1,:))

在这里插入图片描述

使用 odextend 将解扩展到 t f = 35 t_f=35 tf​=35，并将结果添加到原始图中。

sol_new = odextend(sol,@vdp1,35);
x = linspace(20,35,350);
y = deval(sol_new,x);
hold on
plot(x,y(1,:),'r')

在这里插入图片描述

1.3 输入参数

1.3.1 options – options 结构体

options 结构体，指定为结构体数组。使用 odeset 函数创建或修改 options 结构体。有关与每个求解器兼容的选项列表，请参阅 ODE 选项摘要。

示例：
options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot)

指定 1e-5 的相对误差容限、打开求解器统计信息的显示，并指定输出函数 @odeplot 在计算时绘制解。

数据类型： struct

From: ode45-MathWorks

版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请联系我们举报，一经查实，本站将立刻删除。

发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/203735.html原文链接：https://javaforall.net

赞 (0)
0 0

关于作者

全栈程序员-站长的头像

全栈程序员-站长

133.5K 文章
3 粉丝
本网站汇聚当前互联网主流语音，持续更新，欢迎关注公众号“全栈程序员社区”
上一篇 2025年8月9日 下午9:43
下一篇 2025年8月9日 下午10:22


相关推荐

  • win10下pytorch-gpu安装以及CUDA详细安装过程

    win10下pytorch-gpu安装以及CUDA详细安装过程

    win10下pytorch-gpu安装以及CUDA详细安装过程win10下pytorch-gpu安装以及CUDA详细安装过程1.Cuda的下载安装及配置首先我们要确定本机是否有独立显卡。在计算机-管理-设备管理器-显示适配器中，查看是否有独立显卡。可以看到本机有一个集成显卡和独立显卡NVIDIAGetForceGTX1050。接下来，测试本机独立显…

    全栈程序员-站长的头像 全栈程序员-站长
    2026年2月21日
    4
  • Android多线程开发之Callback

    Android多线程开发之Callback

    Android多线程开发之Callback多线程中的返回 Callback 回调 interface

    全栈程序员-站长的头像 全栈程序员-站长
    2026年3月20日
    2
  • win10修改视频默认播放器_win10无法更改默认播放器

    win10修改视频默认播放器_win10无法更改默认播放器

    win10修改视频默认播放器_win10无法更改默认播放器你要用到WindowsPowerShell，它是win10系统自带的一个应用，要打开它，就单击开始菜单中的“所有应用”，然后找到WindowsPowerShell的文件夹，右键单击WindowsPowerShell（注意不是WindowsPowerShellISE），然后单击以管理员身份运行，就打开了。输入以下：get-appxpackage*zunevideo*|r

    全栈程序员-站长的头像 全栈程序员-站长
    2026年4月15日
    6
  • Css Filter Alpha 属性详解

    Css Filter Alpha 属性详解

    Css Filter Alpha 属性详解Alpha 设置透明层次 blur 创建高速度移动效果 即模糊效果 Chroma 制作专用颜色透明 DropShadow 创建对象的固定影子 FlipH 创建水平镜像图片 FlipV 创建垂直镜像图片 glow 加光辉在附近对象的边外 gray 把图片灰度化 invert 反色 light 创建光源在对象上 mask 创建透明掩膜在对象上 shadow 创建偏移固定影子 w

    全栈程序员-站长的头像 全栈程序员-站长
    2026年3月16日
    3
  • Ubuntu 搭建opengrok 流程

    Ubuntu 搭建opengrok 流程

    Ubuntu 搭建opengrok 流程opengrok平台搭建流程软件下载链接:https://pan.baidu.com/s/1kCeXNlj2l3FujyMza3rM0w提取码:iniy搭建前的准备电脑系统电脑系统推荐使用ubuntu16，这版系统较为稳定。细节未更新python环境推荐使用python2.7及以上版本，这一版本相对稳定，python安装细节未更新java环境推荐使用JDK1.8及以上版本，具体安装细节未更新通过java-version和javac-version可以查看版本。Ope

    全栈程序员-站长的头像 全栈程序员-站长
    2022年5月27日
    82
  • 汉澳sinox不受openssl心血漏洞影响并分析修复其漏洞代码

    汉澳sinox不受openssl心血漏洞影响并分析修复其漏洞代码

    汉澳sinox不受openssl心血漏洞影响并分析修复其漏洞代码OpenSSL心血(HeartBleed)漏洞是openssl在2014-04-07发布的重大安全漏洞（CVE-2014-0160）这个漏洞使攻击者可以从server内存中读取64KB的数据，甚至获取到加密流量的密钥。用户的名字和password。以及訪问的内容。主要影响版本号OpenSSL1.0.1到OpenSSL1.0.1f以及OpenSSL1.0….

    全栈程序员-站长的头像 全栈程序员-站长
    2022年7月16日
    28

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注

关注全栈程序员社区公众号