Code stuck at line 210 of interp1 while solving ODE with time-dependent terms.

I'm trying to solve a second-order ODE with ode45. It has two time-dependent tems E and eta (both 1x1000) as well as several constant values A, rho, g, L, and I. Given it's second-order I've created a system of first-order odes. I'm having to interpolate E and eta. I created the following function...

function dqdt = MyODE(t,q,E,eta,time,A,rho,g,L,I)
    E_intp = interp1(time,E,t);
    eta_intp = interp1(time,eta,t);
    dq1dt = q(2);
    dq2dt = (((21*rho*A*g)-(504.*E_intp'.*I.*q(1))-(504.*eta_intp'.*I.*q(2)))./(rho*A*L^4));
    dqdt = [dq1dt; dq2dt];
end

It's used in the main code like this...

t_span = [0 600];
q0 = [0;0];

[t,q] = ode45(@(t,q) MyODE(t,q,E,eta,time,A,rho,g,L,I),t_span,q0);

When I run the code it runs forever and therefore I don't get an error message. It gets to line 210 of the inbuilt interp1 function and stays there. What have I got wrong? Is there an infinite loop?

[Jan]

Remember, that Matlab's ODE integrators are designed to solve smooth functions only. A linear interpolation is not smooth. This can reduce the stepsize to tiny values, such that the accumulated rounding errors dominate the final value. Therefore this is not a reliable way to solve the intergal, although it might produce some fair values by accident.

But stopping inside interp1 is not related to this problem. I assume that you have set a breakpoint by accident. Try dbclear all or restart Matlab.

[Joshua Bowes-Smith]

I'm pretty new to matlab, I'm not really sure what you mean by smooth. Why is there an example with interpolation here: https://uk.mathworks.com/help/matlab/ref/ode45.html#bu3l43b if it's not reliable to use with ode45?

Your suggestions didn't stop interp1 stopping at line 210. Would you like me to attach the complete files to take a look? Please help?

[Torsten]

As already noted, we need to run your code. For this, it is necessary to get the data files for eta and E.

[Joshua Bowes-Smith]

Here you go

[Torsten]

You already posted the code, but we need the data arrays E and eta.

If they are somehow confidential, maybe you can generate artificial data that represent the original data good enough to reproduce the failure of the integrator.

[Jan]

@Joshua Bowes-Smith: "Smooth" means differentiable in any point. A linear interpolation has knees. Then the step size controller of the ODE integrator can fail tremendously. If you are lucky, this produces a warning or error, but sometimes a final value is produced, which should not be called "result", because it can be more or less random.

This is taught in the lessons for numerical maths. The correct solution is to stop and restart the integration at the corresponding points.

Unfortunately the docs of ODE45 do contain an example for a linear interpolation inside the function to be integrated.

[Joshua Bowes-Smith]

There aren't separate E and eta data files. They are simply vectors calculated on lines 32 and 34 of the code. They should appear in the workspace if you force stop the code after running.

[Joshua Bowes-Smith]

How do you stop and restart the integration at the corresponding points?

[Torsten]

Works. Didn't I suggest to try ode15s instead of ode45 ?

Are you sure about the curves for eta and E ?

If you generate eta and E from given formula that depend on time, why don't you directly use these formulae to calculate eta and E in MyODE instead of interpolating them from data arrays ?

Or do you intend to use measurement data later on in your project ?

tmin = 0;
tmax = 600;
N = 1000;
dt = (tmax - tmin)/(N-1);
time = tmin:dt:tmax;

%Newton's cooling parameters
Tenv = 295.15;
T0 = 513.15;
%For a beam of D=2mm, L=50mm, rahman_paper
m = 1.65e-04;
c = 1200;
h = 10;
A = 3.14e-04;
k = (h*A)/(m*c);
%Arrhenius equation parameters
eta_inf = 100;
Ea = 175800;
R = 8.314;
%Young temperature dependence parameters
E0 = 2.3e+09;
B0 = 1.25e+07; %matched to Huang paper 
T_0 = 600; %matched to Huang paper

%Formula 1: Temperautre_Time
T = Tenv + (T0 - Tenv)*exp(-k*time);
%Formula 2: Viscosity_Temperature
eta = eta_inf * exp(Ea./(R*T));
%Formula 3: YoungModulus_Temperature
E = E0 - (B0 .* T .* exp(-T_0./T));

%Figures
figure(1)
plot(time,T,'-r')
grid
xlabel('Time [s]')
ylabel('Temperature [K]')

figure(2)
set(gca,'YScale','log')
plot(time,eta,'-b')
grid
xlabel('Time [s]')
ylabel('Viscosity [Pa.s]')

figure(3)
plot(time,E,'-g')
grid
xlabel('Time [s]')
ylabel('Young Modulus [Pa]')

%ODE parameters
%A = A; already set above
rho = 1050;
g = 9.81;
L = 50e-03;
D = 2e-03;
I = 0.25*pi*(D/2)^4;
t_span = [0 50];
q0 = [0;0];

%Solving ODE
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
%[t,q] = ode15s(@(t,q) MyODE(t,q,E,eta,time,A,rho,g,L,I),t_span,q0);
[t,q] = ode15s(@(t,q) MyODE2(t,q,A,rho,g,L,I),t_span,q0,options);

figure (4)
plot(t,q(:,1),'-o')

function dqdt = MyODE(t,q,E,eta,time,A,rho,g,L,I)
E_intp = interp1(time,E,t);
eta_intp = interp1(time,eta,t);
dq1dt = q(2);
dq2dt = (((21*rho*A*g)-(504.*E_intp'.*I.*q(1))-(504.*eta_intp'.*I.*q(2)))./(rho*A*L^4));
dqdt = [dq1dt; dq2dt];
end

function dqdt = MyODE2(t,q,A,rho,g,L,I)
%Newton's cooling parameters
Tenv = 295.15;
T0 = 513.15;
%For a beam of D=2mm, L=50mm, rahman_paper
m = 1.65e-04;
c = 1200;
h = 10;
%A = 3.14e-04;
k = (h*A)/(m*c);
%Arrhenius equation parameters
eta_inf = 100;
Ea = 175800;
R = 8.314;
%Young temperature dependence parameters
E0 = 2.3e+09;
B0 = 1.25e+07; %matched to Huang paper 
T_0 = 600; %matched to Huang paper
%Formula 1: Temperautre_Time
T = Tenv + (T0 - Tenv)*exp(-k*t);
%Formula 2: Viscosity_Temperature
eta = eta_inf * exp(Ea./(R*T));
%Formula 3: YoungModulus_Temperature
E = E0 - (B0 .* T .* exp(-T_0./T));
dq1dt = q(2);
dq2dt = (((21*rho*A*g)-(504.*E.*I.*q(1))-(504.*eta.*I.*q(2)))./(rho*A*L^4));
dqdt = [dq1dt; dq2dt];
end

[Joshua Bowes-Smith]

Hi again, thanks for this. This seems to be working fine interms of solving the ODE. I'm not getting the results I want though. I think I want to change Ea to a lower value like 12000 however then the graph for q becomes oscillatory which isn't right. Help? It does however settle at a value of 133 which is closer than 8e-09, I want a value of magnitude e+04.

Also why is q 1500953x2 when is should just be solving the ODE for 1000 E and 1000 eta?

[Torsten]

Also why is q 1500953x2 when is should just be solving the ODE for 1000 E and 1000 eta?

All MATLAB integrators do not integrate with a fixed stepsize, but adapt their stepsize during integration. Thus although you prescribe 1000 values for E and eta, it might happen that MATLAB needs 10 or 100000 steps to reach the final time. If you want MATLAB to output the results at prescribed output times, you can specify them in the vector tspan (at the moment, tspan is a two-element vector in your code which means that MATLAB will return the solution for q at all time steps required to reach the final time).