Euler method with tabular inputs and multiple ODE's

I am attempting to model some parameters of a lake, namely algae and P contents, and lake volume.

For this, I have 3 seperate equations (for X, P, and V), and also refer to input data in the form of daily rainfall and evaporation readings that are in an excel spreadsheet.

I'm struggling in how to format the Euler loop(/s?) and getting the data to use the combination of euler interpretation for some variables, and reading from the spreadsheet for each timestep for others.

Fair warning I am very much a beginner user so apologies in advance if my code is messy.

%Load BOM Data
[date, t, Rainfall, Evap, I] = readvars("BOM_Data.xlsx");
nt = size(date,1);

% EULER METHOD
dt = 1;  % step size (days)
time = 0:dt:365;  % the range of time (days)
V = zeros(size(time,1));  % allocate the resulting volume of lake
V(1) = 22410;  % the initial volume of lake -FULL

% Temperature Function
T(t) = 23.1 + 5.8*sin((0.5 + (2*t)/365)*pi);

% Light Intensity function
f_I = (I(t)/I_s) .*exp((-I(t)/I_s)+1);

% EULER METHOD
% Initial conditions and setup
dt = 1;  % step size (days)

% The loop to solve the DE
for t=1:nt
    time(t+1) = time(t) + dt;
    dP = P_in - P_out - r_p;
    dX = r_x - r_d - X_out;
    dV = (Rainfall(t) * A_catch * C / 1000) - Overflow(t) - (Evap(t) * A / 1000);    %units = m3

        r_x(t+1) = Mu_max * f_I * P(t)/(K_s + P(t)) * X(t);
        r_d(t+1) = k_d * Theta^(T - 20) * X(t);
        r_p(t+1) = r_x(t) * 1/Y_xp;

        P_in(t+1) = Rainfall(t) *A_catch * C ;  %units = mg
        P_out(t+1) = P(t) * Overflow(t) / V(t);    %units = mg
        X_out(t+1) = X(t) * Overflow(t) / V(t);    %units = mg
  
    P(t+1) = P(t) + dP * dt;    % mass P in mg
    X(t+1) = X(t) + dX * dt;    % mass X in mg
    V(t+1) = V(t) + dV * dt;    % volume lake in m3
end

Answer 1

[Torsten]

The best way to refer to time-dependent data vectors in integration problems is to use "interp1".

In your case above, e.g., define a function RAINFALL as

RAINFALL = @(time) interp1(t,Rainfall,time)

where t and Rainfall are the arrays you read from "BOM_Data.xlsx" and use this function in the Euler loop as

RAINFALL(time(it))

(You should rename the loop index of the Euler loop form t to it since you overwrite the vector t from the Excel file).

Further obvious errors:

% Temperature Function

T = @(t) 23.1 + 5.8*sin((0.5 + (2*t)/365)*pi);

instead of

T(t) = 23.1 + 5.8*sin((0.5 + (2*t)/365)*pi);

% Light Intensity function

Int = @(time) interp1(t,I,time)
f_I = @(time) (Int(time)/I_s) .*exp((-Int(time)/I_s)+1);

instead of

% Light Intensity function

f_I = (I(t)/I_s) .*exp((-I(t)/I_s)+1);

And in the Euler loop, always refer to the actual index of the variables, e.g.

dP(it) = P_in(it) - P_out(it) - r_p(it);

instead of

dP = P_in - P_out - r_p;