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:
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!