Choosing a numerical method for solving a stiff system of ODES

Hello,

I am in the process of constructing a zero dimensional adiabatic parcel model during which I wish to numerically integrate the system forward in time. The system is a closed system of first order ODEs and it is my understanding that this system is stiff which limits the numerical methods that I can utilise. I am aware that there are built in solvers for stiff problems, I do not wish to use these however I will be comparing the output of my numerical solver with them.

It is also my understanding that for a stiff system of ODEs, implicit methods are required such as the BDF and Adams-moulton solvers. Whilst I have found a large selection of literature on solving non-stiff systems I have struggled to find anything that seems applicable to my problem.

My problem is an initial value problem and the initial conditions are calculated in the main script and are then placed in a state vector of the form:

y_init = [z0, P0, T0, S0, wv0, wl0, r0]; 

Another vector, t, for the time range also exists which is defined in the main script:

t_end = ceil(z_top/V); % choose when to stop the solver (in this case after 250m)

t = 0:h:t_end;         % set the time range, h is the step size

y = zeros(numel(y_init),numel(t)); % initialise an empty matrix

y(:,1) = y_init;       % set the first column values to the initial conditions

At present I have assumed the step size, h, is constant, however this does not have to be the case.

I have written a function which describes the system of ODEs, below is a simplified version of this function with some calculations removed (some of these calculations depend on P, T, S ect...).

function PM_dydt = PM_diff_eq(t,y,r_dry,N,kappa,V)

z = y(1), P = y(2), T = y(3), S = y(4), wv = y(5), wl = y(6), r = y(7);

% 1) Define constants

% CONSTANTS ARE DEFINED HERE

% 2) Calculate parameters

% CALCULATIONS GO HERE

% 3) Define ODE system
 
  dzdt = V;
  dPdt = -1.0*((g*P*V)/(Rd*Tv));
  drdt = (G/r)*(S-Seq);
  dwldt = (4.0*pi*(rho_w/rho_ad))*N*(r^2)*drdt ;
  dwvdt = -1.0*dwldt;
  dTdt = -((g*V)/Cp)-(Lv/Cp)*dwvdt;
  dSdt = ((P/(eps*esw))*dwvdt-((S+1.0)*(((eps*Lv)/(Rd*T^2))*dTdt)+((g*V)/(Rd*T));

  % Create a vector showing the time derivatives
  
  PM_dydt = [dzdt; dPdt; dTdt; dwvdt; dwldt; dSdt; drdt];
  
end

You will notice that many of the derivatives depend on knowing values of another derivative, hence the system has to be solved in the order presented.

Can anyone suggest how to numerically integrate this system forward in time without using the built in solvers? Or offer any advice on how to approach this problem?

Kind regards, Liam

Answer 1

[Drishti]

Hi Liam Holbeche-Smith,

To numerically integrate a stiff system of ODEs without using built-in solvers, you can utilize the implicit Euler method, which is a first-order implicit method suitable for stiff problems.

In this method a system of equations are solved at each time step. In MATLAB, you can utilize the 'fsolve' function from the Optimization Toolbox to do the same.

For better understanding, you can refer to the code snippet provided below:

function implicit_euler_solver
    % Define constants and parameters as per requirement
    
    % Set Initial conditions
 
    % Time-stepping loop
    for i = 1:length(t)-1
        y(:,i+1) = implicit_euler_step(t(i), y(:,i), h, r_dry, N, kappa, V);
    end
end

function y_next = implicit_euler_step(t, y_current, h, r_dry, N, kappa, V)
    % Define the function to solve using fsolve
    function F = implicit_eq(y_next)
        y_dot = PM_diff_eq(t + h, y_next, r_dry, N, kappa, V);
        F = y_next - y_current - h * y_dot;
    end

    % Use fsolve to solve the implicit equation
    options = optimoptions('fsolve', 'Display', 'off');
    y_next = fsolve(@implicit_eq, y_current, options);
end

The above code snippet provides an overview on using euler's method.

Refer to the MathWorks Documentation of 'fsolve' function for more information:

https://in.mathworks.com/help/optim/ug/fsolve.html

I hope this helps in getting started.