1. 全栈程序员必看首页

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

全栈程序员-用户IM • • 未分类

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

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

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

【正版授权，激活自己账号】： Jetbrains全家桶Ide使用，1年售后保障，每天仅需1毛

【官方授权 正版激活】： 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

赞 (0)
全栈程序员-用户IM全栈程序员-用户IM
0 0


相关推荐

  • 达梦数据库的函数_达梦数据库连接命令

    达梦数据库的函数_达梦数据库连接命令

    达梦数据库的函数_达梦数据库连接命令SQL工作笔记-达梦数据库关于时间的函数http://blog.itpub.net/69995127/viewspace-2758308/达梦数据库的查询以及函数的使用

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年10月28日
  • 手机安装busybox教程(去除busybox的xbin连接)

    手机安装busybox教程(去除busybox的xbin连接)

    Busybox,大家应该有所耳闻，其中集成了众多linux的命令和工具。有时候在手机上想使用find，grep等命令但是手机不支持。这时可以在手机上安装busybox.使用方法：1.确定手机型号cpu的指令集，cat/proc/cpuinfo如armv7，armv4,armv6,去网上下载二进制busybox文件。下载地址:http:

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年4月12日
  • 实现dubbo服务降级

    实现dubbo服务降级

    实现dubbo服务降级dubbo降级服务使用dubbo在进行服务调用时，可能由于各种原因（服务器宕机/网络超时/并发数太高等），调用中就会出现RpcException，调用失败。服务降级就是指在由于非业务异常导致的服务不可用时（上面举得例子），可以返回默认值，避免异常影响主业务的处理。dubbo服务降级配置mock配置方式dubbo官方文档上使用一个mock配置，实现服务降级。mock只在出现非业务异常(比如超时，网络

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年6月4日
  • 恩智浦被中国收购(光纤放大器调试图解)

    恩智浦被中国收购(光纤放大器调试图解)

    恩智浦半导体（NXPSemiconductorsN.V.）（纳斯达克：NXPI）日前推出业界首款能进行超宽带应用（470至806MHz）的宽带Doherty功率放大器，新产品采用BLF884P和BLF884PS架构的超宽带Doherty参考设计。全新70WDVB-TLDMOS设计采用能在超宽带范围内工作的恩智浦架构（正在申请专利），为广播发射机带来Doherty拓扑结构的效率提升。

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年4月17日
  • 电脑dnf,DNF卡顿如何解决_DNF卡顿如何解决 教你调整电脑参数畅玩游戏_52PKDNF「建议收藏」

    电脑dnf,DNF卡顿如何解决_DNF卡顿如何解决 教你调整电脑参数畅玩游戏_52PKDNF「建议收藏」

    电脑dnf,DNF卡顿如何解决_DNF卡顿如何解决 教你调整电脑参数畅玩游戏_52PKDNF「建议收藏」DNF卡顿怎么解决？相信很多玩家电脑的配置并不差，但是就是玩DNF会卡。今天就在这里教大家一些优化的方法。让你轻松摆脱DNF卡顿带来的困扰。如果是硬件问题的可以换硬件，如果是软件设置问题的可以优化自己的设置。首先要说的由于系统版本和软件版本的问题，每一项设置带来的提升也会不同，需要各位玩家自己摸索，我也会讲一下自己在不同系统测试的感受。我的电脑信息：下面是方法汇总：一、硬件：1.显卡虽然之前一直有…

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年9月11日
  • idea2019激活码永久_在线激活 idea

    idea2019激活码永久_在线激活

    (idea2019激活码永久)本文适用于JetBrains家族所有ide，包括IntelliJidea，phpstorm，webstorm，pycharm，datagrip等。IntelliJ2021最新激活注册码，破解教程可免费永久激活，亲测有效，下面是详细链接哦~https://javaforall.cn/100143.html…

    全栈程序员-用户IM 全栈程序员-用户IM
    2022年3月28日

发表回复

您的电子邮箱地址不会被公开。

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