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

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!

Answer 1

[Alan Stevens]

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

[Ying Wu]

Hi @Alan Stevens, thanks for your detailed reply! And I think I have solve this problem. But I still have a quick question about the function spline. With the given vectors of dY and P, I can use spline to do cubic spline interpolation by the following commend and get the struct pp. As shown, the range of dY is [0, 0.0763], however, for the dY out of this range, there still exists value of P (see figure below for dY<0 ad dY>0.0763). So my questions is how the function spline to generate the value of P at the dY outside the given range. Thanks!

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];
pp = spline(dY,P);
dY_test = [-0.01:0.01:0.08];
P_test = fnval(pp,dY_test);
plot(dY,P,'ko',dY_test,P_test,'k-');

[Alan Stevens]

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]

Thanks again for your suggestion!