How can I solve this system of ODEs?

1 view (last 30 days)
John
John on 19 Aug 2020
Answered: Sam Chak on 19 Aug 2020
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

Accepted Answer

Sam Chak
Sam Chak on 19 Aug 2020
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.
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')
  2 Comments
John
John on 19 Aug 2020
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!

Sign in to comment.

More Answers (1)

Sam Chak
Sam Chak on 19 Aug 2020
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);

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!