MATLAB ode45的三种传参方式

MATLAB ode45的三种传参方式

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

在以下例子中，我们想要通过ode45求解以下微分方程 { d y 1 d t = − b y 1 y 1 + y 2 + 1 G a d y 2 d t 1 G a d y 2 d t d 2 y 2 d t 2 = b y 1 y 1 + y 2 + 1 G a d y 2 d t d y 2 d t − G a d y 2 d t left{ egin{array}{c} frac{dy_1}{dt}=frac{-by_1}{y_1+y_2+frac{1}{Ga}frac{dy_2}{dt}}frac{1}{Ga}frac{dy_2}{dt}\ \ frac{d^2y_2}{dt^2}=frac{by_1}{y_1+y_2+frac{1}{Ga}frac{dy_2}{dt}}frac{dy_2}{dt}-Gafrac{dy_2}{dt}\ end{array} ight. ⎩⎪⎪⎨⎪⎪⎧dtdy1=y1+y2+Ga1dtdy2−by1Ga1dtdy2dt2d2y2=y1+y2+Ga1dtdy2by1dtdy2−Gadtdy2 而使用ode45求解则需要把上式转换为一阶微分方程组，因而我们作以下换元 令 G a ∗ y 3 = d y 2 d t 令Ga*y_3=frac{dy_2}{dt} 令Ga∗y3=dtdy2 因而上述微分方程组可化为 { d y 1 d t = − b y 1 y 1 + y 3 + 1 G a d y 2 d t 1 G a d y 2 d t d y 3 d t = b y 1 y 3 y 1 + y 2 + y 3 − G a ∗ y 3 d y 2 d t = G a ∗ y 3 left{ egin{array}{c} frac{dy_1}{dt}=frac{-by_1}{y_1+y_3+frac{1}{Ga}frac{dy_2}{dt}}frac{1}{Ga}frac{dy_2}{dt}\ \ frac{dy_3}{dt}=frac{by_1y_3}{y_1+y_2+y_3}-Ga*y_3\ \ frac{dy_2}{dt}=Ga*y_3\ end{array} ight. ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧dtdy1=y1+y3+Ga1dtdy2−by1Ga1dtdy2dtdy3=y1+y2+y3by1y3−Ga∗y3dtdy2=Ga∗y3

可以看到上述微分方程组中有两个常数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