Implementing a code from Berkley Madonna into Matlab

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]

@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.

Answer 1

[William Rose]

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]

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