Euler method with tabular inputs and multiple ODE's
3 次查看(过去 30 天)
显示 更早的评论
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
0 个评论
回答(1 个)
Torsten
2022-10-12
编辑:Torsten
2022-10-12
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;
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!