How can I solve this system of ODEs?
4 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
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
Sam Chak
on 19 Aug 2020
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.
More Answers (1)
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);
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!