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

Question

Ying Wu 2021-8-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/894897-how-to-solve-a-2nd-order-ode-system-with-space-dependent-variable-with-ode45

评论： Ying Wu 2021-8-13

采纳的回答： Alan Stevens

在 MATLAB Online 中打开

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:

I use cubic spline interpolation to fit these points and obtain a piecewise function for P(dY): pp=csapi(dY,P).
Like the time-dependent variable in Matlab ode45 help, I define a vector of dY and obtain the corresponding P(dY).
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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2021-8-9

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/894897-how-to-solve-a-2nd-order-ode-system-with-space-dependent-variable-with-ode45#answer_763517

在 MATLAB Online 中打开

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 个评论
显示 1更早的评论隐藏 1更早的评论

Alan Stevens 2021-8-13

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 2021-8-13