开启掘金成长之旅！这是我参与「掘金日新计划 · 2 月更文挑战」的第 2 天，点击查看活动详情
求解常微分方程常用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;endend
• 这是一个常微分方程系统，其中对自变量时间有不止一个导数。
• 这是一个刚性系统，因为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); % z1z2=z(2); % z2dzdt = [z2 ; -G/L*sin(z1);];endend
