How can I solve this system of ODEs?
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
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
0 comentarios
Respuesta aceptada
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
el 19 de Ag. de 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.
Más respuestas (1)
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);
0 comentarios
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!