MATLAB Answers

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

5 views (last 30 days)
Winner Anigbogu
Winner Anigbogu on 9 Oct 2020
Edited: Winner Anigbogu on 10 Oct 2020
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?

  0 Comments

Sign in to comment.

Answers (1)

John D'Errico
John D'Errico on 9 Oct 2020
Edited: John D'Errico on 9 Oct 2020
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.
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.

  1 Comment

Winner Anigbogu
Winner Anigbogu on 10 Oct 2020
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.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by