Solving Coupled ODE's by ode45

I have several ODE's with initial condition. I've given them as example.

dx/dt = dy/dt + r,

r = a*x+b (a,b are parameters)

dy/dt = z for 1st 30 minutes (z = e*f)

after 30 minutes

dy/dt = z - k*f (k,f are parameters)

How can I solve this?

[James Tursa]

What have you done so far? What specific problems are you having with your code? Typically you would call ode45( ) with a time span of 30 minutes, then take that result as the initial conditions for the next time span and call ode45( ) again with the different derivative function.

[Swachwa Dey]

Thanks for the comment. Actually I am a new user. How can I define the function for 1st ode? I am facing trouble in defining the function for differential equation which has another differential term in the equation.

[James Tursa]

Have you looked at the ode45( ) doc? There are examples there for solving systems of equations:

https://www.mathworks.com/help/matlab/ref/ode45.html

Basically you make a vector (in your case a two element vector since you have two variables x and y) and then write your equations based on this vector.

[Swachwa Dey]

I am not getting something which involves another differential term in an ODE like I said in my post. I can solve something like dx/dt = 3*x+4*y, dy/dt=5*x+9*y. But how can I solve these or define function when the equations look dx/dy = 3*(dy/dt) + 4*x, dy/dt = dz/dt + 4*x and dz/dt = 4*x?

Answer 1

[James Tursa]

Start with your example (assuming that dx/dy was supposed to be dx/dt as you had originally):

dx/dt = 3*(dy/dt) + 4*x, dy/dt = dz/dt + 4*x and dz/dt = 4*x?

Then make these state vector definitions, using Y as the state vector for all the variables:

Y(1) = x

Y(2) = y

Y(3) = z

The derivative function would look something like this

function dYdt = myderiv(t,Y)
x = Y(1);
y = Y(2);
z = Y(3);
dzdt = 4*x; % do this one first, doesn't depend on the others
dydt = dzdt + 4*x; % then this one next since we now have dzdt
dxdt = 3*dydt + 4*x; % then do this one last since we now have dydt
dYdt = [dxdt;dydt;dzdt];
end

[Swachwa Dey]

function graca = gracaALspecies(t, C)

%parametric values
k1f = 4.2e4;
k2f = 4.2e4;
k3f = 5.6e4;
k4f = 2.5e-7;
kcg = 1e-3;
kwf = 1.52e-6;
I = 190;
F = 96485;
Z = 3;
V = 0.006;
K1 = 9.6e-6;
K2 = 5.3e-5;
K3 = 2e-6;
K4 = 2.7e-9;
K_w = 1e-14;
K_s = 4.6e-33;
pH = 7;
epsilon = 1.6;

%variables
c1 = C(1);
c2 = C(2);
c3 = C(3);
c4 = C(4);
c5 = C(5);
c6 = C(6);
c7 = C(7);
c8 = C(8);

rAl = I*epsilon/(F*Z*V);
Al_Sat = K_s*K_w^-3*(10^(-3*pH)+K1*10^(-2*pH)+K2*10^(-pH)+K3+K4*10^pH);

%differential equations
dc8dt = rAl - kcg*(c8-Al_Sat);
dc1dt = dc8dt-k1f*(c1-(c2*c6)/K1);
dc2dt = k1f*(c1-(c2*c6)/K1) - k2f*(c2-(c3*c6)/K2);
dc3dt = k2f*(c2-(c3*c6)/K2) - k3f*(c3-(c4*c6)/K3);
dc4dt = k3f*(c3-(c4*c6)/K3) - k4f*(c4-(c5*c6)/K4);
dc5dt = k4f*(c4-(c5*c6)/K4);
dc6dt = k1f*(c1-(c2*c6)/K1) + k2f*(c2-(c3*c6)/K2)...
    + k3f*(c3-(c4*c6)/K3) + k4f*(c4-(c5*c6)/K4)...
    + kwf*(1-(c6*c7)/K_w);
dc7dt = 3*dc8dt + kwf*(1-(c6*c7)/K_w);
graca = [dc1dt; dc2dt; dc3dt; dc4dt; dc5dt; dc6dt; dc7dt; dc8dt];
end %function file ends here

% I created another script file for follwoing code to run the above file
dt = 10;
tspan = [0 5400];
C0 = [0 0 0 0 0 1e-7 1e-7 0]; %initial values
[t, C] = ode45(@gracaALspecies, tspan, C0, dt)

Sir, I've followed your instruction. I am getting this error message.Could you please let me know where I am making the mistake?