利用ode45求解含控制量并且控制量为离散点的动力学方程

1、写出微分方程函数2、求解

function dy=rigid(t,y)
dy=zeros(3,1);
dy(1)=y(2)*y(3);
dy(2)=-y(1)*y(3);
dy(3)=-0.51*y(1)*y(2);
end
%将微分方程写成函数形式，待调用
options=odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[T Y]=ode45(@rigid,[0 12],[0 1 1],options);
plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'.');
%[0 12]是时间跨度[0 1 1]是初始值，值得注意的是求出的Y为矩阵，画图调用时要注意

在实际的动力学方程中会涉及控制量，而微分方程只是计算状态量，因此计算时需要输入控制量的时间函数，但另一个问题我们计算得到的优化控制量大多时候是离散点，这些点可以拟合或插值得到随时间的函数。

function dydt = myode(t,y,ft,f,gt,g)
%UNTITLED2 ft是f所对应的时间节点，gt是g对应的时间节点，首先对控制量进行插值
%   此处显示详细说明
f=interp1(ft,f,t);
g=interp1(gt,g,t);
dydt=-f.*y+g;
end
%y是状态量，f和g为两个控制量，微分方程写为函数形式调用
clc;clear;close;
%% 控制量数据点
ft=linspace(0,5,25);
f=ft.^2-ft-3;
gt=linspace(1,6,25);
g=3*sin(gt-0.25);
%%  计算，注意@函数句柄调用形式
Tspan=[1 5];%时间跨度
Iy=1;%初始值
[T,Y]=ode45(@(t,y) myode(t,y,ft,f,gt,g),Tspan,Iy);
plot(T,Y);