在ode MATLAB求解器中使用中间值

Question

我正在使用一个僵硬的求解器(ode15s)对一个颂歌系统进行时间积分。它在起作用，但我想加快速度。

方程组以状态空间形式给出：

function [dx] = fun(t,x,M,C,K,other_parameters)
    % Mx'' + Cx' + Kx = F(t)
    % BUNCH OF CALCULATIONS
    F = solveP(x,t);

    A = [zeros(n) eye(n) ; -M\K -M\C];
    b = M\F;

    dx = A*x + b
end

这里的诀窍部分是强迫函数F，它是高度非线性的，依赖于x和t参数。它利用x参数求解Poisson型二维方程（用有限体积法）。力F与泊松方程解成正比。

function [F] = solveP(x,t)
    % initialize solution
    Phi = zeros(Ni,Nj);

    % solve iteratively
    % ...
    while (~converged)
         % some calculations

         % iterative solver
         Phi(i,j) = (aE*Phi(i,j+1) + aW*Phi(i,j-1) + aN*Phi(i+1,j) +...
                     aS*Phi(i-1,j) + S(i,j))/aP;
    end


    % calculate F
    F = sum(Phi(:)); % discrete integration over domain
end

用迭代法求解泊松方程需要一个初始条件，我把它设为零(phi=zeros(Ni,Nj))。我想我可以通过提供一个更好的场的初始估计来提高计算速度（一个更好的初始条件会更快地到达所寻求的答案）。我能想到的最优初始条件（除了=0)是ode求解器前一次迭代（最后一步）中获得的字段的值((k)_initial=(k-1))。

PS：我尝试使用持久变量，但这不是一个好的解决方案。ode解算器在时间前进之前分几点计算函数。persistent变量在每次ode调用odefun时保存收敛字段。这并不完全是我想要的，而且随着时间的推移，这实际上提供了错误的答案。

Answer 1

一种可能的修复方法是在一次迭代后将找到的phi解决方案写入基本工作区。首先，您必须在基本工作区中初始化phi:

Phi = zeros(Ni, Nj);

然后，您可以通过使用evalin在solvep函数中调用这个phi。然后，在找到PHI的更新解决方案后，确保使用assignin在基本工作区中再次分配它。solvep将如下所示：

function [F] = solveP(x,t)
    % initialize solution
    Phi = evalin('base','Phi');

    % solve iteratively
    % ...
    while (~converged)
         % some calculations

         % iterative solver
         Phi(i,j) = (aE*Phi(i,j+1) + aW*Phi(i,j-1) + aN*Phi(i+1,j) +...
                     aS*Phi(i-1,j) + S(i,j))/aP;
    end


    % calculate F
    F = sum(Phi(:)); % discrete integration over domain

    % assign updated Phi in base workspace
    assignin('base', 'Phi2', Phi)
end