Effacer les filtres
Effacer les filtres

Solve a system of two differential equations

3 vues (au cours des 30 derniers jours)
Pedro Gouveia
Pedro Gouveia le 11 Août 2020
Modifié(e) : hosein Javan le 11 Août 2020
Hello,
I am trying to solve this set of differential equations in order to time in MATLAB but i have no clue how to do it without using methods like Euler, Runge-Kuta and with those I am not getting valid results. The equations are the following:
Everything but and H and its dependent equations are known fixed values. The goal is to determine how these two change in order to time. One thing is that it is important to know the instant when .
I am still new to MATLAB so i have no clue how to use a numerical solver like ode45. I have made research here in the forum but still can't figure it out. I would appreciate all the help. Thanks in advance!

Réponses (1)

hosein Javan
hosein Javan le 11 Août 2020
Modifié(e) : hosein Javan le 11 Août 2020
generally in order to solve a system of ode s; you must first define your equation in the form of (Xdot = f(X,U)).
"U" is the input function that must be defined and is time dependant. in your case everything is well-defined except input function. however I managed to do a code which would probably help. this code does not work correctly until you define parameters. the function "Vw(H)" is not defined also.
% define constants:
Aout = 1;
Vb = 1;
Vw0 = 1;
r = 1;
rho_w =1;
P1 = 1;
Patm = 1;
a = 1;
g = 1;
% define functions:
Vw = @(x) x;
A = @(X) arrayfun (@(x) 0.177*x + 0.007, X);
B = @(X) arrayfun (@(H) integral(@(Z)Aout./A(Z),0,H), X);
% define ode
% assume: x(1) = uout, x(2) = H :
F = @(t,x) [
(-0.5*((Aout/A(x(2)))^2-1)*x(1)^2 - ((Vb-Vw0)/(Vb-Vw(x(2))))^r*P1/rho_w + Patm/rho_w - (a+g)*x(2))/B(x(2))
Aout*x(1)/A(x(2))
];
% ode options:
tss = 1e-4; % time stepsize in seconds
tend = 0.02; % end time of simulation in seconds
tspan = 0:tss:tend;
x0 = [0,0]; % initial condition
options = odeset('AbsTol',1e-3,'MaxStep',tss);
% solve ode:
[t,x] = ode23tb(F, tspan, x0, options);
x(:,1) = uout;
x(:,2) = H;
plot(t,uout,t,H);

Community Treasure Hunt

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

Start Hunting!

Translated by