MATLAB ode45的三种传参方式

MATLAB ode45的三种传参方式

使用ode45函数时，需要自定义函数。而有的时候自定义函数会有一些内参，如何传参便成为了问题。

在以下例子中，我们想要通过ode45求解以下微分方程
$$\{\begin{array}{c}\frac{d{y}_1}{dt}=\frac{-b{y}_1}{y_1+y_2+\frac{1}{Ga}\frac{d{y}_2}{dt}}\frac{1}{Ga}\frac{d{y}_2}{dt}\\ \\ \frac{d^2{y}_2}{d{t}^2}=\frac{b{y}_1}{y_1+y_2+\frac{1}{Ga}\frac{d{y}_2}{dt}}\frac{d{y}_2}{dt}-Ga\frac{d{y}_2}{dt}\end{array}$$
而使用ode45求解则需要把上式转换为一阶微分方程组，因而我们作以下换元
$$令Ga\ast y_3=\frac{d{y}_2}{dt}$$
因而上述微分方程组可化为
$$\{\begin{array}{c}\frac{d{y}_1}{dt}=\frac{-b{y}_1}{y_1+y_3+\frac{1}{Ga}\frac{d{y}_2}{dt}}\frac{1}{Ga}\frac{d{y}_2}{dt}\\ \\ \frac{d{y}_3}{dt}=\frac{b{y}_1{y}_3}{y_1+y_2+y_3}-Ga\ast y_3\\ \\ \frac{d{y}_2}{dt}=Ga\ast y_3\end{array}$$

可以看到上述微分方程组中有两个常数b和Ga，因而我们想求解不同的b和Ga下的y值，但是传统的ode45中规定我们方程传入参数只能有t和y，那我们便可以采用以下三种方式把常量作为参数传入其中

方法一

首先便是可以直接利用ode45给的传参输入变量来输入传参
即ode45(函数句柄，起止时间（两个数）/指定时间（向量），初始值，flag（置空即可），传参)

c=jet(5);
b=0.1;
hold off
for i = 1:5
    ga=0.01*i;
    [t,y]=ode45(@df1,[0 300],[1000 1 0],[],b,ga);
    plot(t,y(:,2),'Color',c(i,:) ...
        ,'DisplayName',['ga=',num2str(0.1*i)] ...
        ,'LineWidth',2)
    hold on
end
legend show

function dy=df1(t,y,b,ga)
    dy=zeros(3,1);
    dy(1)=-b*y(1)*y(2)/sum(y);
    dy(2)=b*y(1)*y(2)/sum(y)-ga*y(2);
    dy(3)=ga*y(2);
end

方法二（推荐）

其次便是可以利用函数句柄进行传参
这种传参方式可以扩展到其他需要输入自定义函数的MATLAB函数（比如fmincon，其规定自定义函数传入参数只能有y），而且比起方法三来说，这种方法不会因为使用全局变量而产生一些难以发现的错误

c=jet(5);
b=0.1;
hold off
for i = 1:5
    ga=0.01*i;
    [t,y]=ode45(@(t,y) df1(t,y,b,ga),[0 300],[1000 1 0]);
    plot(t,y(:,2),'Color',c(i,:) ...
        ,'DisplayName',['ga=',num2str(0.1*i)] ...
        ,'LineWidth',2)
    hold on
end
legend show

function dy=df1(t,y,b,ga)
    dy=zeros(3,1);
    dy(1)=-b*y(1)*y(2)/sum(y);
    dy(2)=b*y(1)*y(2)/sum(y)-ga*y(2);
    dy(3)=ga*y(2);
end

方法三

最后便是可以利用全局变量进行传参
注意：本方法不能在实时脚本中使用

c=jet(5);
global b ga
b=0.1;
hold off
for i = 1:5
    ga=0.01*i;
    [t,y]=ode45(@df2,[0 300],[1000 1 0]);
    plot(t,y(:,2),'Color',c(i,:) ...
        ,'DisplayName',['ga=',num2str(0.1*i)] ...
        ,'LineWidth',2)
    hold on
end
legend show

function dy=df2(t,y)
    dy=zeros(3,1);
    dy(1)=-b*y(1)*y(2)/sum(y);
    dy(2)=b*y(1)*y(2)/sum(y)-ga*y(2);
    dy(3)=ga*y(2);
end