How to solve a 2nd-order ODE system with space-dependent variable with ODE45?

3 vues (au cours des 30 derniers jours)
Ying Wu
Ying Wu le 9 Août 2021
Commenté : Ying Wu le 13 Août 2021
Hi, I use ode45 function to solve the following 2nd-order ODE system:
where X(t) and Y(t) are function of time, and the parameter P is dependent on the first dirivative of Y(t).
I have 15 scattered points for P and dY:
dY = [0, 0.0026, 0.0082, 0.0138, 0.0194, 0.0250, 0.0306, 0.0362, 0.0418, 0.0475, 0.0532, 0.0589, 0.0647, 0.0705, 0.0763]
P = [0, 0.0423, 0.1534, 0.2336, 0.2915, 0.3628, 0.4166, 0.4587, 0.5572, 0.8002, 0.8204, 0.5559, 0.3184, 0.1405, -0.0603]
My main steps include:
  1. I use cubic spline interpolation to fit these points and obtain a piecewise function for P(dY): pp=csapi(dY,P).
  2. Like the time-dependent variable in Matlab ode45 help, I define a vector of dY and obtain the corresponding P(dY).
  3. Treat (dY,P) as inpus of odefun and use interp1 to update P with dY at each time step during the numerical computation.
The key code are shown below:
% cubic spline interpolation
csi = csapi(dY,P);
dY = [0:0.0001:0.1]';
P = fnval(csi,dy);
% ode45
[time,sol] = ode45(@(t,Y) odefun(t,Y,dY,P),[0 1000],[0,0,0.01,0);
function dYdt = odefun(t,Y,dY,P)
CCFy = interp1(ddY,P,abs(Y(4)));
% Y(1) = X
% Y(2) = X'
% Y(3) = Y
% Y(4) = Y'
dYdt = zeros(4,1);
dYdt(1) = Y(2);
dYdt(2) = 0.51*(0.03*Y(1) + 0.0025*P - 0.045*Y(4) - Y(3)) + 0.887*Y(4) + ...
1.1642*(1-1175.51*Y(1)^2)*Y(2)- 1.774*Y(1);
dYdt(3) = Y(4);
dYdt(4) = 0.03*Y(1) + 0.0025*P - 0.045*Y(4) - Y(3);
end
But the results are incorrect. Could anyone please give me some help? Thanks!

Réponse acceptée

Alan Stevens
Alan Stevens le 9 Août 2021
More like this perhaps
dY = [0, 0.0026, 0.0082, 0.0138, 0.0194, 0.0250, 0.0306, 0.0362, 0.0418, 0.0475, 0.0532, 0.0589, 0.0647, 0.0705, 0.0763];
P = [0, 0.0423, 0.1534, 0.2336, 0.2915, 0.3628, 0.4166, 0.4587, 0.5572, 0.8002, 0.8204, 0.5559, 0.3184, 0.1405, -0.0603];
tspan = [0 1000];
ic = [0,0.01,0,0];
[time,sol] = ode45(@(t,Y) odefun(t,Y,dY,P),tspan,ic);
subplot(2,1,1)
plot(time,sol(:,1)),grid
xlabel('time'),ylabel('Y')
subplot(2,1,2)
plot(time,sol(:,3)),grid
xlabel('time'),ylabel('X')
function dYdt = odefun(~,Y,dY,P)
% Y(1) = Y
% Y(2) = Y'
% Y(3) = X
% Y(4) = X'
p = @(dydt) spline(dY,P,dydt);
dYdt = zeros(4,1);
dYdt(1) = Y(2);
dYdt(2) = 0.03*Y(3) + 0.0025*p(Y(2)) -0.045*Y(2);
dYdt(3) = Y(4);
dYdt(4) = 0.51*dYdt(2) + 0.887*Y(2) - 1.774*Y(3) + ...
1.1642*(1-1175.51*Y(3)^2)*Y(4);
end
  3 commentaires
Alan Stevens
Alan Stevens le 13 Août 2021
My advice is NEVER extrapolate outside the range of the known data unless you have a physical model for the behaviour (which the spline fit doesn't give you)!
Ying Wu
Ying Wu le 13 Août 2021
Thanks again for your suggestion!

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Produits


Version

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by