# Implementing a code from Berkley Madonna into Matlab

2 views (last 30 days)
Justin Lee on 8 Apr 2021
Commented: Cris LaPierre on 8 Apr 2021
I am trying to implement a code from Berkley Madonna into Matlab. I want to carry out a simulation and produce the results graphically by using ode solvers directly. Here is this Berkley Madonna code I am trying to incorporate:
METHOD RK4
STARTTIME = 0
STOPTIME = 300 {minutes}
DT = 0.02
Km = 0.0184
Vg = 2.4
Vl = 1.08
Vc = 11.56
Vm = 25.76
FL = 1.35
Fm = 0.95
D=0.2
{D is a function of 5g ethanol dose and 100 kg body weight}
ks = -0.049*(D) + 0.0545
INIT Vs = 1
d/dt(Vs) = -ks*Vs
INIT Cg = 0
INIT Cl = 0
INIT Cc = 0
INIT Cm = 0
Vmax = 0.1012
d/dt(Cg) =((((2/3)*FL)*(Cc - Cg) +(ks*Vs))/Vg)
rel = (Vmax*Cl)/(Km + Cl)
d/dt(Cl) = ((FL*(0.33*Cc + .66*Cg - Cl) - rel)/Vl)
d/dt(Cc) = ((-FL*(Cc - Cl) - Fm*(Cc - Cm))/Vc)
d/dt(Cm) = (Fm*(Cc-Cm)/Vm)
William Rose on 8 Apr 2021
@Justin Lee, I'm not familiar with Berkley Madonna and I bet most other readers aren;t. If you post the equations you are trying to simulate, instead of a bunch of code that people cannot parse, you will be more likely to get a helpful reply.
You can format your equations by clicking the capital sigma icon about the message window. The equation editing window which pops up uses Latex formatting. I always thought Latex was too complicated for me, but I finally gave it a try a week ago, and it is easier than I expected. Every time I want to do something, I do a google search: "Latex for subscript", "Latex for a fraction", etc. Example (I think this is one of your differential equations): I recommend reading Intro to solving ODEs in Matlab and the Matlab help for ode45(). One of the things you will learn is that you don't need DT=0.02 (which your code defines) because ode45() is a 4th order Runge Kutta routine with adaptive stepsize - the routine figures out and adjusts the step size as it goes.

William Rose on 8 Apr 2021
Edited: William Rose on 8 Apr 2021
Justin,
The attached code shows a way to solve differential equaitons like yours in Matlab, and plot the results. The output is shown below. You should check that the diff eqs in
function dxdt=JustinsODE(t,x)
are what you want, and that all constants are correct.
>> JustinLeeDiffEq
>> Cris LaPierre on 8 Apr 2021
Since I spent time on this, I figured I'd share my solution as well. I have converted Berkley Madonna simulations to MATLAB simulations before, and the references pointed to are great, as BM is just a differential equation solver.
I formatted the plot to appear similar to what is generated in BM
STARTTIME = 0;
STOPTIME = 300;
% initial values [Vs Cg Cl Cc Cm]
y0 = [1;0;0;0;0];
% solve ode
[t,y] = ode45(@odefun,[STARTTIME STOPTIME],y0);
% Plot results
yyaxis left
plot(t,y(:,[1,3]))
ylabel("Vs,Cl")
yyaxis right
plot(t,y(:,[2,4]))
ylabel("Cg,Cc")
legend(["Vs","Cg","Cl","Cc"]) function ddt = odefun(t,y)
% Constants
Km = 0.0184;
Vg = 2.4;
Vl = 1.08;
Vc = 11.56;
Vm = 25.76;
FL = 1.35;
Fm = 0.95;
D=0.2;
ks = -0.049*D + 0.0545;
Vmax = 0.1012;
% Current concentrations
Vs = y(1);
Cg = y(2);
Cl = y(3);
Cc = y(4);
Cm = y(5);
% Differential equations
ddt(1,1) = -ks*Vs; % Vs
ddt(2,1) =((((2/3)*FL)*(Cc - Cg) +(ks*Vs))/Vg); % Cg
rel = (Vmax*Cl)/(Km + Cl);
ddt(3,1) = ((FL*(0.33*Cc + .66*Cg - Cl) - rel)/Vl); % Cl
ddt(4,1) = ((-FL*(Cc - Cl) - Fm*(Cc - Cm))/Vc); % Cc
ddt(5,1) = (Fm*(Cc-Cm)/Vm); % Cm
end