我正在使用刚性求解器(ode15s)对ODE系统进行时间积分。它工作正常,但我想加快速度。
方程组以状态空间形式给出:
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参数来求解泊松型二维方程(通过有限体积法)。力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))的前一次迭代(最后一步)中获得的φ场的值。
底线是:如何在ode解决方案中使用/保存中间值?
PS:我尝试使用持久变量,但这不是一个好的解决方案。 ode求解器在推进之前计算几个点的函数。每次ode调用odefun乐趣时,持久变量都会保存收敛的φ字段。这不完全是我想要的,随着时间的推移,这实际上提供了错误的答案。