ode45函数
ode45实际上是数值分析中数值求解微分方程组的一种方法,4阶五级Runge-Kutta算法。
调用方法
其实这种方程的每一个状态变量都是t 的函数,我们可以从现代控制理论的状态空间来想。因此返回[ t , x ],其中t是一个列向量,x是n × n的矩阵,它的每一列就是其中一个状态变量随t的变化值。
Fun就是你要求解的微分方程组,微分方程组必须化成现代控制理论中的一阶微分方程组形式,之后定义。m函数亦或是匿名函数f = @ ( t , x )来描述微分方程组的右边部分,并且是列向量。(注意匿名函数一定是先t后x并且一定有t占位)
tspan就是你要求解的t区间,x 0 就是初始状态。options是一些选项。
例子
Rossler微分方程组
f = @(t,x)[-x(2)-x(3);x(1)+0.2*x(2);0.2+(x(1)-5.7)*x(3)];
[t,y] = ode45(f,[0,100],[0;0;0]);
plot(t,y);
figure;
plot3(y(:,1),y(:,2),y(:,3))
Riccati 微分方程
P1 = [1 0 0;0 3 0;0 0 5];
[t,p] = ode45(@riccati,[0.5,0],P1(:));
plot(t,p)
function U=riccati(t,x)
A = [6 6 17;1 0 -1;-1 0 0];
B = [0 0 0;0 4 2;0 2 1];
C = [1 2 0;2 8 0;0 0 4];
p = [x(1),x(2),x(3);x(4),x(5),x(6);x(7),x(8),x(9)];
K = A'*p+p*A+p*B*p+C;
U=[K(1,:)';K(2,:)';K(3,:)'];
end
注意这个地方给的是t=0.5s时的状态值,但是Runge-Kutta算法可以倒着算,因此tspan写成[0.5,0],这个时候的初始值就是P1(0.5)。
设置ode45精度
主要是绝对精度和相对精度,设置方式如下:
clc;clear;close all;
f = @(t,x)[x(1)*(x(1)^2+x(2)^2-2)-4*x(1)*x(2)^2;4*x(1)^2*x(2)+x(2)*(x(1)^2+x(2)^2-2)];
option = odeset('RelTol',1e-8,'AbsTol',[1e-8;1e-8]);
[t,y] = ode45(f,[0,10],[0.5;0.5],option);
plot(t,y(:,1),t,y(:,2))
RelTol为相对精度,一维数据,AbsTol相对精度,要给符合状态变量的维数。之后把option传入ode45即可。