I'm using ode15s and i keep getting a Warning regarding integration tolerance

I'm using ode15s and i keep getting a Warning regarding integration tolerance. I'm not sure what I'm doing wrong or what I need to change...

This is my set of code:

--------------------------------------------------------
% File Name: ROM_Fringe3.m

N = 3; % Terms (Modes of Vibration)

w1 = 3.51562; % 1st Natural Frequency of Cantilever Beam
c1 = -0.734; % 1st Constant Coefficient

w2 =  22.0336; % 2nd Natural Frequency of Cantilever Beam
c2 = -1.0185; % 2nd Constant Coefficient

w3 = 61.70102; % 3rd Natural Frequency of Cantilever Beam
c3 = -0.992; % 3rd Constant Coefficient

b = 0.001; % Damping Coeficient
d = 0.1; % Voltage
f = 0.26; % Fringe Effect

sigma_pt3 = -0.00168908; % Frequency Test Point (-0.00109406 -0.00117026 -0.00128908 -0.00148374 -0.00168908 -0.00204088 -0.00231234)
ohm = sigma_pt3 + w1; % NOTE: Ohm = Omega + Sigma

phi1 = @(z) -1.*(cos(sqrt(w1).*z) - cosh(sqrt(w1).*z) + c1.*(sin(sqrt(w1).*z) - sinh(sqrt(w1).*z))); % Cantilever Mode Shape 1
phi2 = @(z) -1.*(cos(sqrt(w2).*z) - cosh(sqrt(w2).*z) + c2.*(sin(sqrt(w2).*z) - sinh(sqrt(w2).*z))); % Cantilever Mode Shape 2
phi3 = @(z) -1.*(cos(sqrt(w3).*z) - cosh(sqrt(w3).*z) + c3.*(sin(sqrt(w3).*z) - sinh(sqrt(w3).*z))); % Cantilever Mode Shape 3

tspan = 0:0.05:15000;
y0 = [0;0.3;0;0.3;0;0.3]; % [Init. Amplitude; Init. Velocity] N = 1:[0.2;0]

options = odeset('Mass', @MassMatrix3_Fringe); %  Mass Matrix Option (File: MassMatrix3_Fringe)
[t, y] = ode15s(@ODEFunc3_Fringe,tspan,y0,options,ohm); %Try ode23t, ode15s, or ode45

u1 = y(:,1);
u2 = y(:,3);
u3 = y(:,5);

u = u1.*phi1(1) + u2.*phi2(1) + u3.*phi3(1);

plot(t,u,'b','LineWidth',0.25);
hold on; grid on;
yline(0,'k','LineWidth',2) % Base Line
xlabel('Time, t')
ylabel('Amplitude, u')
[top, bottom] = title("ROM Time Response", "b = " + b + "  f = " + f + "  \delta = " + d + "  \omega_1 = " + w1 + "  \sigma = " + sigma_pt3 + "  N = " + N + "");
bottom.FontAngle = 'italic';
ylim([-2 2]);
legend('u','Base Line')

--------------------------------------------------------

This is the error I'm getting:

>> ROM_Fringe3

Warning: Failure at t=4.370445e+03. Unable to meet integration tolerances without reducing the step size below the

smallest value allowed (1.455192e-11) at time t.

> In ode15s (line 655)

In ROM_Fringe3 (line 38)

Answer 1

[Walter Roberson]

The code in ODEFunc3_Fringe encounters a singularity in the system at about 4370.4 seconds.

That code (which you do not show) might have a bug.

The code in MassMatrix3_Fringe (which you do not show) might have a bug.

Your initial conditions might be wrong.

The system might be inherently unstable.

As this appears to be homework: sometimes the point of telling you to "Try ode23t, ode15s, or ode45" is to show you that sometimes it matters which function you call, that sometimes calling one will lead to a singularity when the others might not. So it is not certain that you are doing anything "wrong": this just might happen to be the expected outcome with that ode*() routine.

[Miguel Martinez]

When I try ode45, I get this:

Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.

> In odemassexplicit>ExplicitSolverHandleMass34 (line 67)

In odefcncleanup>@(t,y)oldFcn(t,y,inputArgs{:}) (line 11)

In ode45 (line 311)

In ROM_Fringe3 (line 38)

When I try ode23t, it continues to run on and on for over an hour and never plots anything.

[Miguel Martinez]

The code for ODEFunc3_Fringe:

function ODE3Fringe = ODEFunc3_Fringe(t,y,ohm)

format compact; 
% File Name: ODEFunc3_Fringe.m

N = 3; % Terms (Modes of Vibration)

b = 0.001; % Damping Coefficient
d = 0.1; % Voltage
f = 0.26; % Fringe Effect

w1 = 3.51562; % 1st Natural Frequency of Cantilever Beam
w2 =  22.0336; % 2nd Natural Frequency of Cantilever Beam
w3 = 61.70102; % 3rd Natural Frequency of Cantilever Beam

% G-Coefficients
g1 = 0.782976200483;
g11 = 0.999986079475;
g111 = 1.477828935063;
g1111 = 2.348985866193;

g2 = 0.433618702528;
g12 = -0.000498667205;
g121 = -0.399376993670;
g21 = -0.000498667205;
g22 = 1.000319912250;
g211 = -0.399376993670;
g212 = 0.955061090714;
g221 = 0.955061090714;
g222 = 0.452886076782;

g3 = 0.433618702528;
g13 = -0.000498667205;
g23 = 1.000319912250;
g31 = -0.000498667205;
g32 = 1.000319912250;
g33 = 1.000319912250;
g131 = -0.399376993670;
g231 = 0.955061090714;
g232 = 0.452886076782;
g311 = -0.399376993670;
g321 = 0.955061090714;
g331 = 0.955061090714;
g312 = 0.955061090714;
g322 = 0.452886076782;
g332 = 0.452886076782;
g313 = 0.955061090714;
g323 = 0.452886076782;
g333 = 0.452886076782;

B1 = -b.*y(2).*g11 + b.*y(2).*y(1).*g111 - (w1.^2).*y(1).*g11 + (w1.^2).*y(1).*y(1).*g111 + f.*d.*((cos(ohm.*t)).^2).*g1; 
B2 = -b.*(y(2).*g21 + y(4).*g22) + b.*(y(2).*y(1).*g211 + y(4).*y(3).*g222) - ((w1.^2).*y(1).*g21 + (w2.^2).*y(3).*g22) + ((w1.^2).*y(1).*y(1).*g211 + (w2.^2).*y(3).*y(3).*g222) + f.*d.*((cos(ohm.*t)).^2).*g2; 
B3 = -b.*(y(2).*g31 + y(4).*g32 + y(6).*g33) + b.*(y(2).*y(1).*g311 + y(4).*y(3).*g322 + y(6).*y(5).*g333) - ((w1.^2).*y(1).*g31 + (w2.^2).*y(3).*g32 + (w3.^2).*y(5).*g33) + ((w1.^2).*y(1).*y(1).*g311 + (w2.^2).*y(3).*y(3).*g322 + (w3.^2).*y(5).*y(5).*g333) + f.*d.*((cos(ohm.*t)).^2).*g3; 
ODE3Fringe = [y(2); B1; y(4); B2; y(6); B3];
end

[Miguel Martinez]

the code for MassMatrix3_Fringe:

function MassMx3_F = MassMatrix3_Fringe(t,y,ohm)

format compact; 
% File Name: MassMatrix3_Fringe.m

N = 3; % Terms (Modes of Vibration)

% G-Coefficients
g1 = 0.782976200483;
g11 = 0.999986079475;
g111 = 1.477828935063;
g1111 = 2.348985866193;

g21 = -0.000498667205;
g22 = 1.000319912250;
g12 = -0.000498667205;
g121 = -0.399376993670;
g211 = -0.399376993670;
g212 = 0.955061090714;
g221 = 0.955061090714;
g222 = 0.452886076782;

g13 = -0.000498667205;
g23 = 1.000319912250;
g31 = -0.000498667205;
g32 = 1.000319912250;
g33 = 1.000319912250;
g131 = -0.399376993670;
g231 = 0.955061090714;
g232 = 0.452886076782;
g311 = -0.399376993670;
g321 = 0.955061090714;
g331 = 0.955061090714;
g312 = 0.955061090714;
g322 = 0.452886076782;
g332 = 0.452886076782;
g313 = 0.955061090714;
g323 = 0.452886076782;
g333 = 0.452886076782;


A11 = g11 - y(1).*g111;
A12 = g12 - y(1).*g121;
A13 = g13 - y(1).*g131;

A21 = g21 - (y(1).*g211 + y(3).*g212);
A22 = g22 - (y(1).*g221 + y(3).*g222);
A23 = g23 - (y(1).*g231 + y(3).*g232);

A31 = g31 - (y(1).*g311 + y(3).*g312 + y(5).*g313);
A32 = g32 - (y(1).*g321 + y(3).*g322 + y(5).*g323);
A33 = g33 - (y(1).*g331 + y(3).*g332 + y(5).*g333);

MassMx3_F = zeros(6,6);
MassMx3_F(1,1) = 1;
MassMx3_F(2,2) = A11;
MassMx3_F(2,4) = A12;
MassMx3_F(2,6) = A13;
MassMx3_F(3,3) = 1;
MassMx3_F(4,2) = A21;
MassMx3_F(4,4) = A22;
MassMx3_F(4,6) = A23;
MassMx3_F(5,5) = 1;
MassMx3_F(6,2) = A31;
MassMx3_F(6,4) = A32;
MassMx3_F(6,6) = A33;

end

[Torsten]

Do the results look reasonable at t=4370 ?

[Miguel Martinez]

They do, yes.

The program plots up until that time, t, then stops.

[Walter Roberson]

ODE3Fringe = [y(2); B1; y(4); B2; y(6); B3];

That second last position, the one that shows as y(6) here:

If you record time and those values when time starts exceeding 4370, then at least for ode15s, as you get close to the singularity, the absolute values fluctuate around +/- 2e-7 ... but sometimes (but not consistently!) go out to about +/- 4E-3 near the singularity. That is a high relative change compared to the change in inputs, so ode15s decides the equation has become unstable.

The second last is just being copied from the last input, so it is not immediately obvious what the problem might be.

If you plot the values of the last output near the singularity, nothing odd shows up -- though it does appear to be accelerating (or I am reading the plot wrong.)

[Walter Roberson]

I made a change to the objective function. Near the top I added

if t > 4300; y = sym(y, 'd'); end

and after ODE3Fringe is assigned to with the current code, I added

ODE3Fringe = double(ODE3Fringe);

Notice this does not change the precision of input or output, only the precision used for the calculations within any one call.

This in addition to the code I had already added to record outcomes near the singularity:

    if t >= 4370
        global mem
        mem(end+1,:) = [t; y; ODE3Fringe];
    end

which I had before the double()

The resulting code took much longer because of the symbolic operations -- and it did nearly 1000 more iterations before it gave up at t=4.372225e+03 (instead of t=4.370445e+03)

The output plots for the recorded intervals look significantly different than before, with several outputs showing a more cyclic appearance.

The third output, y(4) in your ODE3Fringe variable, shows a near singularity near the original problem time, but with the additional resolution it manages to power through that. But at the new problem time, there is a second singularity, and this time that output gets fairly steep, but along a curve. At the same time, the fifth output, y(6) in your ODE3Fringe variable, suddenly gets quite steep, but more like a jump.

My speculation is that you are being limited by precision -- but that there is also something about the derivatives that is driving to effective singularity anyhow. I predict that if higher precision were carried (somehow) that maybe you would get through the new oscillation, but that the next one (probably at time related to 2*pi/rho later) would be even worse.