How can I solve this system of ODEs?

1 visualización (últimos 30 días)
John
John el 19 de Ag. de 2020
Respondida: Sam Chak el 19 de Ag. de 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

Respuesta aceptada

Sam Chak
Sam Chak el 19 de Ag. de 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 comentarios
Sam Chak
Sam Chak el 19 de Ag. de 2020
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
John el 19 de Ag. de 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!

Iniciar sesión para comentar.

Más respuestas (1)

Sam Chak
Sam Chak el 19 de Ag. de 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);

Categorías

Más información sobre Programming en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by