matlab中通过ode函数求解常微分方程附加简单的钟摆模型
matlab中通过ode函数求解常微分方程附加简单的钟摆模型
求解常微分方程常用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