Unable to meet integration tolerances without reducing the step size below the smallest value allowed

I'm trying to solve an N+n coupled equation simultaneusly, where N is the total number of Mode shapes and n is the number of resonators I have on the Euler Bernoulli beam. The two equations look pretty straight forward, but I've been unable to get any result as it is showing me the above error. I believe I'm getting something wrong or probably need to add something to the ODEsolver. Does anyone have an idea where the problem is from, please?

Below is my function file, which I created separately:

function dpdt = myfunction3(t,p,w,prm)
Fo = prm.Fo;
% Zi = prm.Zi;
mr = prm.mr;
nBt = prm.nBt;
xrM = prm.xrM;
BtM2 = prm.BtM2;
Kv = prm.Kc;
Kd = prm.Kd;
n = prm.nxr;
N = prm.nBt;
mysum1 = zeros(1,nBt);
mysum2 = zeros(1,nBt);
mysum3 = zeros(1,n);
mysum4 = zeros(1,n);
acrtn1 = nan(1,nBt);

for j = 1:N
for i = 1:n
    mysum1(j) = mysum1(j) + (Phi(xrM(i),BtM2(j),prm).*Kv.*(p(1+2*i)-Phi(xrM(i),BtM2(j),prm).*p(1).*p(1))); %kv is very stiff
    mysum2(j) = mysum2(j) + (Phi(xrM(i),BtM2(j),prm).*Kd.*(p(2+2*i)-Phi(xrM(i),BtM2(j),prm).*p(2)));
end

acrtn1(j) = -(Cj(BtM2(j),prm).*p(2)+Kj(BtM2(j),prm).*p(1)-...
    (Phi(0,BtM2(j),prm).*Fo.*cos(w.*t)+mysum1(j)+mysum2(j)))./(Mj(BtM2(j),prm)); %First equation, Mj is the mass, Kj stiffness of beam
%Cj is 2*etha*sqrt(Mj,Kj)
end

acrtn2 = nan(1,n);
for i = 1:n 
for j = 1:N
   mysum3(i) = Phi(xrM(i),BtM2(j),prm).*p(1);
   mysum4(i) = Phi(xrM(i),BtM2(j),prm).*p(2);
end
acrtn2(i) = -(Kd.*(p(2+2*i)-mysum4(i)) + Kv.*(p(1+2*i)-mysum3(i)))./mr; %Second equation
end


dpdt = nan(2*(N+n),1);
for j = 1:N
    [dpdt(2*j-1),dpdt(2*j)] = deal(p(2*j),acrtn1(j));
end
for i = 1:n
    [dpdt(2*N+2*i-1),dpdt(2*N+2*i)] = deal(p(2*N+2*i),acrtn2(i)); 
end

end

This is my odesolver file:

%%
clearvars
clc
close all
tic
count = 0;
parm2  % I created a struct file for my parameters
xM = [0,prm.L];
nx = length(xM);
%%%% Mode
BtM = [4.73/l,7.853/l,10.995/l,14.137/l,17.2787/l,20.39/l,23.5/l,26.63/l,29.81/l];
BtM2= BtM/prm.nxr;
prm.BtM2 = BtM2;
nBt = length(BtM2);
prm.nBt = nBt;
%%%% Frequency
wM = 2*pi.*linspace(1,1000,1000); %Change back the frequency to 1000
nw = length(wM);
prm.wM = wM;
%%%% prepare for loops
mystep = 2*pi/(20*max(wM));
t2 = (0:mystep:1)';
num = length(t2);

%% Find pM
pM = nan(nw,num,2*(nBt+nxr));

for j = 1:nw
    [t,p] = ode45(@(t,p) myfunction3(t,p,wM(j),prm),[0 1],zeros(1,2*(nBt+nxr)));
    pM(j,:,:) = interp1(t,p,t2);
    count = count + 1;
    prcnt = count/(nw);
end
toc
save('genDispDat.mat','pM','-v7.3')
save('genDispDat.mat','t2','wM','num','nx','nw','-append')

%% Find W(x,t)
load('genDispDat.mat','pM')
Wxt = nan(nw,nx,num);
for j = 1:nw
    mymat = zeros(nx,num);
    for i = 1:nBt
        mymat = mymat + (Phi(xM,BtM2(i),prm)')*squeeze(pM(j,:,2*i-1));
    end
    Wxt(j,:,:) = mymat;
end
disp('Done W(x,t), Saving Wxt data...')
save('genDispDat.mat','Wxt','-append')

Plot function:

x0 = 1;
Bt = 2;
freq = 8;
close all
set(gcf,'unit','normalize','position',[0,0,1,1],'color','w')
load('genDispDat.mat','pM','Wxt','myamp','fresp','t2','wM','trans_s')

if x0 == 1
    mytext_x = '0';
elseif x0 == 2
    mytext_x = 'L';
else
    error('the first argument can only be 1 or 2')
end
subplot(2,2,1)
hold on
plot(t2,squeeze(pM(freq,:,1)),'-','linewidth',2)
plot(t2,squeeze(pM(freq,:,2)),'--','linewidth',2) %just interested in seeing a few plots of P(t)
legend('1','2','3','4')
xlabel('t(s)')
ylabel('p(t)')
title({'p(t)',['\beta = ',num2str(Bt),', f = ',num2str(freq),' Hz']})
set(gca,'FontWeight','Bold','FontSize',15)
grid on

subplot(2,2,2)
plot(t2,squeeze(Wxt(freq,1,:)));
hold on
plot(t2,squeeze(Wxt(freq,2,:)));

Can anyone pickout something wrong with this, please?

Answer 1

[John D'Errico]

Anytime you see this error message:

"Unable to meet integration tolerances without reducing the step size below the smallest value allowed"

This suggests you need to use a stiff solver. ODE45 is not the correct code to be using in those cases. Use ODE15s or ODE23s as good choices. They are designed to handle stiff sytems.

https://en.wikipedia.org/wiki/Stiff_equation

The idea is ODE45 is an adaptive solver. When it hits a difficult point in the solve, it reduces the step size. When it can no longer reduce the step size below a point on the order of eps, it must then give up, as it would take infinitely long to solve the numerical problem as posed. Then it reports the error you saw.

[Winner Anigbogu]

I actually tried it using ODE15s too, it gave me the same reply. I think the challenge is with my coupling of the N+n equations. I somehow believe I'm not getting something right there. I'm trying to solve J(1,2,3,..,N) and i(1,2,3,...,n) equations simultaneously.