Use intermediate value on ode MATLAB solver

Question

I am doing a time integration of a system of ODEs using a stiff solver (ode15s). It is working, but I want to speed things up.

The system of equations is given in state space form:

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

The trick part here is the forcing function F. It is highly nonlinear and depends on the x and t parameters. It uses the x parameters to solve a Poisson-type 2D equation (by the Finite Volume method). The force F is proportional to the Poisson equation solution.

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

Solving the Poisson equation by a iterative method requires an initial condition, which I set to zero (Phi=zeros(Ni,Nj)). I thought I could improve the speed of calculations by providing a better initial estimative of the ϕ field (a better initial condition would take faster to reach the sought answer). The optimal initial condition I can think (besides ϕ=0) is the value of the ϕ field obtained in the previous iteration (the last step) of the ode solver (ϕ(k)_initial=ϕ(k-1)).

Bottom line is: how do I use/save intermediate values in the ode solution?

PS: I tried using the persistent variables, but that is not a good solution. The ode solver calculates the function in several points before advancing in time. The persistent variable saves the converged ϕ field every time the ode calls the odefun fun. That is not exactly what I want and this actually provides wrong answers as time marching goes on.

From what you show of solveP, F is only dependent on the number of iterations (and thus not on x). If that is the case, why not first determine Phi for the required number of iterations, then the cumsum gives F for every time step. So this would just be a lookup table. Am I missing something? — rinkert
Are you looking for the method to save F variable, the same as you save dx? Do I understand it correctly? — Karls
@rinkert actually F is dependent on x. Changed the code to show that. The code is more extensive than that. The Poisson-type equation is a diffusion equation. The dependency of the equation with the x variable appears both in the diffusion coefficient and the source term. I solve the diffusion equation by a finite voluem method (in this particular case, the final equations are equal to a finite difference discretization of the partial differential equation). — Thales
@Karls yes. As the time marches, the differences of subsequente Phi fields are not too big. Since I solve the discretized form of the pde by an iterative method, if I could provide a better initial estimative for the Phi field (instead of initializing with zeros in every ode iteration), I believe I could save some time. If I were to write an ode solver, I could pass is it as input and output arguments - but I will never be able to write an ode solver nearly as good as MATLAB's solvers. — Thales

rinkert rinkert · Accepted Answer · 2019-04-08T16:11:53

A possible fix is to write the found solution for Phi to the base workspace after one iteration. First you have to initialize Phi in the base workspace:

Phi = zeros(Ni, Nj);

Then you can just recall this Phi in the solveP function by using evalin. Then after you found an updated solution for Phi, be sure to assign it again in the base workspace using assignin. solveP would then look like this:

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

Use intermediate value on ode MATLAB solver

1 Answers