matlab中通过ode函数求解常微分方程附加简单的钟摆模型

✅作者简介：热爱科研的算法开发者，Python、Matlab项目可交流、沟通、学习。

????个人主页：算法工程师的学习日志

求解常微分方程常用matlab中的ode函数，该函数采用数值方法用于求解难以获得精确解的初值问题。ODE是一个包含一个独立变量（例如时间）的方程以及关于该自变量的一个或多个导数。在时域中，ODE是初始值问题，因此所有条件在初始时间t=0指定。

Matlab有几个不同的函数（内置）用于ODEs的解决方案。这些解算器可以与以下语法一起使用：

[outputs] = function_handle(inputs)
[t,state] = solver(@dstate,tspan,ICs,options)

其中state-一个数组，ODE的解（每次状态的值）。

solver-求解器函数，比如ode45、ode23等

dstate- 包含求导公式的函数句柄

tspan- 时间范围，比如[0,5]

ICs- 求解变量的初始状态

options-其他配置参数，比如rtol、atol等

积分器使用我们已经知道并重复的信息计算y(t)的附近值。

高阶数值方法以速度为代价减少误差：

•欧拉方法-一阶展开

•中点法-二阶扩展

•Runge Kutta-四阶扩展

几种不同的求解器对比

[t,state] = ode45(@dstate,tspan,ICs,options)计算步骤：

1.在一个文件中定义tspan、IC和选项（例如call_dstate.m） ，用来设置ode45

2.在另一个文件中定义常量和求导数（例如dstate.m）或作为调用内的函数dstate

3.运行call_dstate.m

4.将结果进行分析

举个例子：

function [t,y] = call_dstate()
tspan = [0 9]; % 时间范围
y0 = 10; % 初始值
% 调用ode
[t,y] = ode45( @dstate ,tspan ,y0);
plot(t,y)
disp([t,y]) % 输出结果
function dydt = dstate (t,y)
    alpha=2; gamma=0.0001;
    dydt = alpha* y-gamma *y^2;
end
end

• 这是一个常微分方程系统，其中对自变量时间有不止一个导数。

• 这是一个刚性系统，因为y1和y2变化剧烈，因此我们需要ode15。

•这次我们将为调用函数（call_osc.m）和ode函数（osc.m）创建单独的文件

为了模拟这个系统，创建一个包含方程的函数osc。

方法1：在列向量中预先分配空间，并填充导数函数

function dydt = osc(t,y)
 dydt = zeros(2,1)
 dydt(1) = y(2);
 dydt(2) = 1000*(1 - y(1)^2)*y(2) - y(1);
end

方法2：对微分函数进行矢量化

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

现在用上述初始条件在0到3000的时间间隔内求解。创建y1随时间的散点图。

function [T,Y] = call_osc()
 tspan = [0 3000];
 y1_0 = 2;
 y2_0 = 0;
 [T,Y] = ode15s(@osc,tspan,[y1_0 y2_0]);
 plot(T,Y(:,1),'o')
end

对于一个简单的钟摆模型​

它的数学模型为：

令：

，则

function [] = call_pend()
tspan=[0 2*pi]; % 
z0=[pi/3,0]; % 初始数值
[t,z]=ode23(@pend,tspan,z0);
plot(t,z(:,1))
function dzdt = pend(t,z)
G=9.8; L=2; % 常数
z1=z(1); % z1
z2=z(2); % z2
dzdt = [z2 ; -G/L*sin(z1);];
end
end