# Choosing a numerical method for solving a stiff system of ODES

1 view (last 30 days)
Liam Holbeche-Smith on 20 Jan 2021
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);
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