1 view (last 30 days)

Show older comments

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!

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

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

Start Hunting!