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

时间:2022-07-04 00:54:30

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

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

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

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

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等

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

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

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

•欧拉方法-一阶展开

•中点法-二阶扩展

•Runge Kutta-四阶扩展

几种不同的求解器对比



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

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

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

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

3.运行call_dstate.m

4.将结果进行分析

举个例子:

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

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

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

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

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

• 这是一个刚性系统,因为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

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

对于一个简单的钟摆模型

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

它的数学模型为:

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

令:

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

,则

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

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

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