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

lichengwan

求解常微分方程常用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])% 输出结果functiondydt = dstate (t,y)alpha=2; gamma=0.0001;dydt = alpha* y-gamma *y^2;endend




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

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

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

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

方法1：在列向量中预先分配空间，并填充导数函数
functiondydt = osc(t,y)dydt = zeros(2,1)dydt(1) = y(2); dydt(2) = 1000*(1 - y(1)^2)*y(2) - y(1);end
方法2：对微分函数进行矢量化
functiondydt = 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); % z1z2=z(2); % z2dzdt = [z2 ; -G/L*sin(z1);];endend