MATLAB Answers

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

1 view (last 30 days)
Ying Wu
Ying Wu on 9 Aug 2021
Commented: Ying Wu on 13 Aug 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!

Accepted Answer

Alan Stevens
Alan Stevens on 9 Aug 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 Comments

Sign in to comment.

More Answers (0)

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by