How can I solve this system of ODEs?

I am trying to program a model in MatLab, which consists of a system of three ODEs. I've put in all of the inputs, equations and tried to use the "syms" and "dsolve" codes to deal with the system of ODEs, but I can't seem to get what I want from the model. Specifically, I would ideally like to solve for "O" (second to last equation) over model time, given initial values of O = 6 and I = 8.5. My issues (from what I can tell) are twofold:

1) How can I specify these initial conditions of O and I?

2) How can I write the code for a function that gives me t and O?

Thank you!

% Inputs %

t = linspace(0,1000,1);      % model time 
V1 = 2.489*10^16;            % vol NA (km^3)
V2 = 3.75*10^15;             % vol NSe (km^3)
V3 = 9.018*10^14;            % vol AO (km^3) 
F2 = 0.072;                  % Fw rate NS (Sv)
F3 = 0.064;                  % Fw rate AO (Sv)
S1 = 35.2;                   % NA S
S2 = 34.9;                   % NS S
S3 = 34;                     % AO S
kO = 10714;                  % Hydraulic constant O (Sv)
kE = 3472;                   % Hydraulic constant E (Sv)
a = 10^-4;                   % Thermal expansion coefficient (per K)
B = 8*10^-4;                 % Haline contraction coefficient (per psu)
DT = 8;                      % T DIFF NA - Others

% SYSTEM OF DIFFERENTIAL EQUATIONS %

syms S1(t) S2(t) S3(t);

ode1 = diff(S1) == -1/2*((E+abs(O)+abs(I))*(S1-S2))-(E*(S2-S3))+F2+F3;
ode2 = diff(S2) == 1/2*((E+abs(O)+abs(I))*(S1-S2))-F2;
ode3 = diff(S3) == (E*(S2-S3))-F3;
odes = [ode1;ode2;ode3];
    
S = dsolve(odes);
S1sol(t) = S.S1;
S2sol(t) = S.S2;
S3sol(t) = S.S3;

% Specify Initial Conditions %
cond1 = S1(0) == 35.2;
cond2 = S2(0) == 34.9;
cond3 = S3(0) == 34;
conds = [cond1;cond2;cond3];
[S1sol(t), S2sol(t), S3sol(t)] = dsolve(odes, conds);

% Vol * dS/dt %
ans1 = V1*S1sol(t);
ans2 = V2*S2sol(t);
ans3 = V3*S3sol(t);

% Calculations %

E = kE*(B*(S2-S3));          % Vol Transport E
O = kO*(a*DT-B*(S1-S2));     % VOL TRANSPORT O
I = E+O;                     % Vol Transport I

Answer 1

[Sam Chak]

Try this, but please check if the ODEs are entered correctly. You should use the method described in the link. I used this method because I don't want to save the odefcn m-file in my folder.

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

clear all; clc
tspan  = 0:0.01:20;
y0 = [35.2; 34.9; 34];
	V1 = 2.489*10^16;	% vol NA (km^3)
	V2 = 3.75*10^15;	% vol NSe (km^3)
	V3 = 9.018*10^14;	% vol AO (km^3) 
	F2 = 0.072;		% Fw rate NS (Sv)
	F3 = 0.064;		% Fw rate AO (Sv)
	kO = 10714;		% Hydraulic constant O (Sv)
	kE = 3472;		% Hydraulic constant E (Sv)
	a = 10^-4;		% Thermal expansion coefficient (per K)
	B = 8*10^-4;		% Haline contraction coefficient (per psu)
	DT = 8;			% T DIFF NA - Others
[t, y] = ode45(@(t,y) [-1/2*(((kE*(B*(y(2) - y(3)))) + abs(kO*(a*DT - B*(y(1) - y(2)))) + abs(kE*(B*(y(2) - y(3))) + kO*(a*DT - B*(y(1) - y(2)))))*(y(1) - y(2))) - ((kE*(B*(y(2) - y(3))))*(y(2) - y(3))) + F2 + F3; 
		1/2*(((kE*(B*(y(2) - y(3)))) + abs(kO*(a*DT - B*(y(1) - y(2)))) + abs(kE*(B*(y(2) - y(3))) + kO*(a*DT - B*(y(1) - y(2)))))*(y(1) - y(2))) - F2; 
		((kE*(B*(y(2) - y(3))))*(y(2) - y(3))) - F3], tspan, y0);

figure(1)
plot(t, y)
grid on
xlabel('Time, s')
ylabel('States, x(t)')
legend('S1', 'S2', 'S3')
title('Time responses of S1, S2, S3')

kO = 10714;				% Hydraulic constant O (Sv)
kE = 3472;				% Hydraulic constant E (Sv)
a = 10^-4;				% Thermal expansion coefficient (per K)
B = 8*10^-4;				% Haline contraction coefficient (per psu)
DT = 8;					% T DIFF NA - Others
E = kE*(B*(y(:, 2) - y(:, 3)));		% Vol Transport E
O = kO*(a*DT-B*(y(:, 1) - y(:, 2)));	% Vol Transport O
I = E + O;				% Vol Transport I

figure(2)
plot(t, E)
grid on
xlabel('Time, s')
ylabel('E')
title('Time response of Vol Transport E')

figure(3)
plot(t, O)
grid on
xlabel('Time, s')
ylabel('O')
title('Time response of Vol Transport O')

figure(4)
plot(t, I)
grid on
xlabel('Time, s')
ylabel('I')
title('Time response of Vol Transport I')

[Sam Chak]

Hi John

are the functions of the states {}:

,

,

.

In my codes, the entire set of ODEs is expressed in terms of {}. Thus, if you solve the ODEs for , , , then you can find , , . Note that the initial conditions of depend on the initial conditions of . You can check the plots above.

[John]

Hi Sam,

Thank you very much. I noticed that and deleted my former comment. Really great solution to this headache!

This is definitely working as it should now, but I am missing one part of the left side of each equation. Rather than the left side being:

dy(1)/dt

dy(2)/dt

dy(3)/dt

It should be a constant (V) multiplied by the derivatives as:

V1*dy(1)/dt

V2*dy(2)/dt

V3*dy(3)/dt

Thus, I am wondering if it is possible to either build this constant into the model, or extract vectors from from ode45 solutions and multiply them by these constants?

Thanks again for your help!

Answer 2

[Sam Chak]

Thank you for your acceptance. V1, V2, V3 are just constants. You can move V1, V2, V3 to the right-hand side of the ODEs by doing the division. The amended code should look like this:

[t, y] = ode45(@(t,y) [(-1/2*(((kE*(B*(y(2) - y(3)))) + abs(kO*(a*DT - B*(y(1) - y(2)))) + abs(kE*(B*(y(2) - y(3))) + kO*(a*DT - B*(y(1) - y(2)))))*(y(1) - y(2))) - ((kE*(B*(y(2) - y(3))))*(y(2) - y(3))) + F2 + F3)/V1); 
		(1/2*(((kE*(B*(y(2) - y(3)))) + abs(kO*(a*DT - B*(y(1) - y(2)))) + abs(kE*(B*(y(2) - y(3))) + kO*(a*DT - B*(y(1) - y(2)))))*(y(1) - y(2))) - F2)/V2; 
		(((kE*(B*(y(2) - y(3))))*(y(2) - y(3))) - F3)/V3], tspan, y0);