How to solve a 2nd-order ODE system with space-dependent variable with ODE45?
1 次查看(过去 30 天)
显示 更早的评论
Hi, I use ode45 function to solve the following 2nd-order ODE system:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/706822/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/706827/image.png)
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 个评论
采纳的回答
Alan Stevens
2021-8-9
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 个评论
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)!
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!