加快matlab for循环

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

1. 接下来，注意不必在 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 方法。

• 其次，在不知道您要查找的事件的情况下，也许您可​​以使用对数搜索而不是线性搜索。例如二分法。这将大大减少您必须解方程的次数。