加速matlab循环

Question

我有一个包含非线性项的5个ODE系统。我试图在一些范围内改变3个参数，以查看哪些参数将产生我正在寻找的必要行为。 问题是我用3个for循环编写了代码，并且需要很长时间才能获得输出。 当我遇到满足ODE事件的参数集时，我也将参数值存储在循环中。

这就是我在matlab中实现它的方法。

function [m,cVal,x,y]=parameters()

b=5000;
q=0;
r=10^4; 
s=0;
n=10^-8;

time=3000;

m=[];
cVal=[];
x=[];
y=[];

val1=0.1:0.01:5;
val2=0.1:0.2:8;
val3=10^-13:10^-14:10^-11;


for i=1:length(val1)
    for j=1:length(val2)
        for k=1:length(val3)
        options = odeset('AbsTol',1e-15,'RelTol',1e-13,'Events',@eventfunction);
        [t,y,te,ye]=ode45(@(t,y)systemFunc(t,y,[val1(i),val2(j),val3(k)]),0:time,[b,q,s,r,n],options);

         if length(te)==1
             m=[m;val1(i)];
             cVal=[cVal;val2(j)];
             x=[x;val3(k)];
             y=[y;ye(1)];

         end   
        end

    end
end

有没有其他方法可以用来加快这个过程？

个人资料查看器结果

我用简单的格式编写了ODE系统

function s=systemFunc(t,y,p)
        s= zeros(2,1);
        s(1)=f*y(1)*(1-(y(1)/k))-p(1)*y(2)*y(1)/(p(2)*y(2)+y(1));
        s(2)=p(3)*y(1)-d*y(2);
end

f，d，k是常数参数。 方程式比这里的方法更复杂，因为它是一个由5个ODE组成的系统，其中许多非线性项相互作用。

Answer 1

托马索是对的。预分配将节省一些时间。

但我猜想，由于你在循环中运行ode45，所以从根本上没有太多可以做的事情。 ode45本身可能是瓶颈。

我建议你分析你的代码，看看瓶颈在哪里：

profile on 
parameters(... )
profile viewer 

我猜想ode45是个问题。可能你会发现你应该把时间集中在优化systemFunc代码的性能上。但是在运行探查器之前你不会知道。

编辑

根据分析器输出和附加代码，我看到一些有用的东西

1. 看起来你的价值观的矢量化正在伤害你。代替 @（T，Y）systemFunc（T，Y，[VAL1（i）中，val2的（j）中，VAL3（K）]）

尝试

@(t,y)systemFunc(t,y,val1(i),val2(j),val3(k))

您的系统功能定义为

function s=systemFunc(t,y,p1,p2,p3)
        s= zeros(2,1);
        s(1)=f*y(1)*(1-(y(1)/k))-p1*y(2)*y(1)/(p2*y(2)+y(1));
        s(2)=p3*y(1)-d*y(2);
end

2. 接下来，请注意您不必在systemFunc中预先分配空间，只需将它们组合在输出中：

    function s=systemFunc(t,y,p1,p2,p3)
       s = [ f*y(1)*(1-(y(1)/k))-p1*y(2)*y(1)/(p2*y(2)+y(1)), 
             p3*y(1)-d*y(2) ];
    end

最后，请注意ode45内部约占运行时间的1/3。没有太多你能做到的。如果你能忍受它，我建议你将'AbsTol'和'RelTol'增加到更合理的数字。这些值非常小，并且正在使ode45运行很长时间。如果你可以忍受它，尝试将它们增加到1e-6或1e-8之类的东西，看看性能提高了多少。或者，根据功能的平滑程度，您可以使用不同的积分器（如ode23）做得更好。但您的里程将根据问题的顺利程度而有所不同。

Answer 2

我有两个建议给你。

1. 预分配存储结果的向量，并使用增加的索引将它们填充到每次迭代中。
2. 由于您使用的选项始终相同，因此仅在循环外实例化一次。

最终代码：

function [m,cVal,x,y] = parameters()

    b = 5000;
    q = 0;
    r = 10^4; 
    s = 0;
    n = 10^-8;
    time = 3000;
    options = odeset('AbsTol',1e-15,'RelTol',1e-13,'Events',@eventfunction);

    val1 = 0.1:0.01:5;
    val1_len = numel(val1);
    val2 = 0.1:0.2:8;
    val2_len = numel(val2);
    val3 = 10^-13:10^-14:10^-11;
    val3_len = numel(val3);
    total_len = val1_len * val2_len * val3_len;

    m = NaN(total_len,1);
    cVal = NaN(total_len,1);
    x = NaN(total_len,1);
    y = NaN(total_len,1);
    res_offset = 1;

    for i = 1:val1_len
        for j = 1:val2_len
            for k = 1:val3_len
                [t,y,te,ye] = ode45(@(t,y)systemFunc(t,y,[val1(i),val2(j),val3(k)]),0:time,[b,q,s,r,n],options);

                if (length(te) == 1)
                    m(res_offset) = val1(i);
                    cVal(res_offset) = val2(j);
                    x(res_offset) = val3(k);
                    y(res_offset) = ye(1);
                end

                res_offset = res_offset + 1;
            end
        end
    end

end

如果您只想保留已正确计算的结果值，则可以删除函数底部包含NaNs的行。对其中一个向量进行索引就足以清除所有内容：

rows_ok = ~isnan(y);
m = m(rows_ok);
cVal = cVal(rows_ok);
x = x(rows_ok);
y = y(rows_ok);

Answer 3

继续其他建议，我还有2条建议：

• 您可能想尝试使用不同的解算器，ODE45用于非僵硬问题，但从外观来看，您的问题可能看起来很僵硬（参数具有不同的数量级）。尝试使用ode23s方法。
• 其次，在不知道您正在寻找哪个事件的情况下，也许您可​​以使用对数搜索而不是线性搜索。例如二分法。这将严重削减你必须解决方程的次数。