加速 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

我有两个建议给你。

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

最终代码：

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 2

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

但我猜想从根本上说，您无能为力，因为您处于运行ning ode45 循环中。 ode45 本身可能就是瓶颈。

我建议您分析代码以查看瓶颈所在：

profile on 
parameters(... )
profile viewer

我猜 ode45 是问题所在。您可能会发现实际上应该将时间集中在优化 systemFunc 代码以提高性能上。但是在运行分析器之前，您不会知道这一点。

编辑

根据探查器输出和附加代码，我发现了一些有用的东西

您的值的向量化似乎正在伤害您。而不是

@(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

接下来，注意你不必在systemFunc中预先分配space，只需在输出中组合它们：

 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 3

在其他建议的延续下，我还有 2 个建议给您：

您可能想尝试使用不同的求解器，ODE45 适用于非刚性问题，但从外观上看，您的问题可能是刚性的（参数顺序不同量级）。例如，尝试使用 ode23s 方法。
其次，在不知道您要查找哪个事件的情况下，也许您可以使用对数搜索而不是线性搜索。例如二分法。这将大大减少您必须求解方程的次数。

加速 matlab for 循环

Speeding up matlab for loop

performance

matlab

for-loop

differential-equations