Warning: Failure at t=1.130644e-05. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.710505e-20) at time t.
2 views (last 30 days)
Show older comments
ANDREAS CHARALAMBOUS
on 14 Aug 2018
Commented: ANDREAS CHARALAMBOUS
on 14 Aug 2018
Hello everyone,i am trying to simulated a PSA process with matlab.I am a really new in matlab.I start from trying to simluated the adsorption step which is my first step,i get thiw error and i dont really interstand why.Any help will be much appreciated.I attach the model which trying to simulating
if true
%PROCESS FOR AIR SEPARATION.THE MODEL IS IN DIMENSIONLESS FORM.
% NUMBER OF STEPS:
% ADSORPTION:1
% CO-CURRENT DEPRESSURIZATION:2
% COUNTE CURRENT BLOWDOWN:3
% PURGE:4
% PRESSUREIZATION:5
%for nstep=1:5 %if nstep==1 %ADSORPTION STEP xl=0; xr=1; n=61; dx=(xr-xl)/(n-1); velo_f=0.23; %m/s; P=1.;%bar; Tf=298;%K R=8.314;%J/mol*K P=1.;%BAR
T0=Tf; Tw=Tf;
th0=T0/Tf; thw= Tw/Tf;
%Mixture composition at inlet and at start %inlet x=0 yf1=0.78; yf2=1-yf1;
%at start t=0 y10=0; y20=1-y10; P10=P*y10; P20=P*y20;
%Air:78% N2(1),22% O2(2) cf=1.e5*(P/R*Tf); %mol/m^3;
c1f=cf*yf1/cf; c2f=cf*yf2/cf;
c10=cf*y10/cf; c20=cf*y20/cf;
% Mixture: 1=N2,2=O2 MW1=0.028; %kg/mol MW2=0.032; %kg/mol
%Mixture transport properties in the bulk at P,T conditions DM0=1.5e-5; %m2/s Pref=1.5; %ref.pressure DM=DM0*Pref/P; %m2/s,scales with 1/P mu=1.846e-5; % Pa*s,assumed constant with P
%Fixed bed properties eb=0.4 ; L=2.5; %m Rbed=0.5; % m
%Pellet structural properties ep=0.39; Rp=0.475e-3; % m rho_particle =720/(1-eb); %kg/m3 so that rho_bulk =(1-eb)*rho_particle=720 kg/m3
cp_gas=995.3406; %heat capacity of air J/kgK cp_solid=1171.575; %heat capacity of the solid J/kgK lamda=0.418; %solid thermal conductivity
area_bed=2/Rbed; hw=0;
%---------------------------------------------------------------------------------------- % Equilibrium adsorption isotherm data for air
bp1a=(1.05e-3)/1.e5; % attention bp1a must be in 1/Pa in matlab code
bp1b=2012.92;
bp2a=(3.51e-3)/1.e5;
bp2b=1544.43;
qm1a=1.4796;
qm1b=-0.00167;
qm2a=0.3182;
qm2b=-0.000145;
% Enthalpies of adsorption
DH1=-18033.89029; %J/mol not cal/mol DH2=-13222.06341; %J/mol not cal/mol
% Determine Langmuir parameters at T=T0 bp10=bp1a*exp(bp1b/T0) ; bp20=bp2a*exp(bp2b/T0) ;
qm10=qm1a+qm1b*T0; qm20=qm2a+qm2b*T0;
Ke10=qm10*bp10 ; Ke20=qm20*bp20 ;
% Determine Langmuir parameters at T=Tf bp1f=bp1a*exp(bp1b/Tf) ; bp2f=bp2a*exp(bp2b/Tf) ;
qm1f=qm1a+qm1b*Tf; qm2f=qm2a+qm2b*Tf; Ke1f=qm1f*bp1f ; Ke2f=qm2f*bp2f ;
AMW=yf1*MW1+yf2*MW2; %kg/mol rho_gas=cf*AMW; tr=L/velo_f;
qe1f = Ke1f*(P*10^5)*yf1 /( 1 + bp1f*(P*10^5)*yf1 + bp2f*(P*10^5)*yf2 ); qe2f = Ke2f*(P*10^5)*yf2 /( 1 + bp1f*(P*10^5)*yf1 + bp2f*(P*10^5)*yf2 );
qe10 = Ke10*(P*10^5)*y10 /( 1 + bp10*(P*10^5)*y10 + bp20*(P*10^5)*y20 ); qe20 = Ke20*(P*10^5)*y20 /( 1 + bp10*(P*10^5)*y10 + bp20*(P*10^5)*y20 ); %---------------------------------------------------------------------------- De1=2.22e-6; De2=De1*sqrt(MW1/MW2);
Re=(rho_gas*(eb*velo_f)*2*Rp)/mu; Sc=mu/((rho_gas*DM)); Sh=2+1.1*(Sc^(1/3))*(Re^(0.6)); kf=(Sh*DM)/(2*Rp); tdif1=1/(15*ep*De1/Rp^2); tdif2=1/(15*ep*De2/Rp^2); tfilm=Rp/(3*kf); k1inv=tdif1+tfilm; k1=1/k1inv; k2inv=tdif2+tfilm; k2=1/k2inv;
%k1=1000; %k2=1000;
%"LINE DRIVING FORCE RATE APPROXIMATION GAS SEPARATION BY ADSORPTION %PROCESSES",RALPH T.YANG,pg.296,eq.8.47 kc1=15*De1/(Rp^2); kc2=15*De2/(Rp^2); %---------------------------------------------------------------------------- %kc1 = 1000; % N2 mass transfer in the micropores (zeolite crystals) %kc2 = 1000; % O2 mass transfer in the micropores (zeolite crystals)
% Check these correlations
Dz = velo_f/667; lamz= Dz*cp_gas; lame=lamz+lamda/eb;
%tdisp = (Dz/uf^2)*(1-eb)/eb
% Here we non-dimensionalize wrt tr=L/u0 Pe=velo_f*L/Dz; Peth=rho_gas*cp_gas*velo_f*L/lame;
% I put these values to comply with Rege and Fortran file Pe=150; Peth=2; %-------------------------------------
Bi=hw*area_bed*tr/eb*rho_gas*cp_gas; Le=((1-eb)/eb)*rho_particle*cp_solid/(rho_gas*cp_gas); F=((1-eb)/eb)*(tr/c1f)*(3/Rp)*(R/P); A1=(1/(1+Le))+((1-eb)/eb)*(1/(rho_gas*cp_gas*Tf))*qe1f; A2=Bi/(1+Le); A3=(1/(1+Le))+((1-eb)/eb)*(1/(rho_gas*cp_gas*Tf))*qe20; dPdt=0;
param =[Pe Peth Dz DM A1 A2 F rho_particle rho_gas k1 k2 kc1 kc2 ep tr... DH1 DH2 eb Le n thw L cp_gas Rp c1f Bi dPdt bp1a bp2a bp1b bp2b... Tf qm1a qm2a qm1b qm2b qe1f qe2f qe10 qe20 A3 velo_f P cf R]';
for i=1:n c10(i)=0.78; velo_0(i)=1; cp10(i)=c10(i); cp20(i)=1-c10(i); q10(i) =qe10/qe1f; q20(i) =qe20/qe20; th0(i) =T0/Tf; end
u0=[c10 velo_0 cp10 cp20 q10 q20 th0].';
for i=1:n if i==1 x(i) = 0; else x(i) = x(i-1)+dx ; end end t0=0; %lower bound for time tf=10; %upper bound for time m=101; %time points tspan=linspace(t0,tf,m); [t,dudt]=ode15s(@(t,u)ADSORPTION_fun(u,dx,param ),tspan,u0); function [ dudt ] = ADSORPTION_fun(u,dx,param ) %UNTITLED7 Summary of this function goes here %----------------------------------------------------------- %PARAMETERS Pe= param(1,1); Peth=param(2,1); Dz =param(3,1); DM =param(4,1); A1 =param(5,1); A2 =param(6,1); F =param(7,1); rho_particle =param(8,1); rho_gas=param(9,1); k1=param(10,1); k2=param(11,1); kc1=param(12,1); kc2=param(13,1); ep=param(14,1); tr=param(15,1); DH1=param(16,1); DH2=param(17,1); eb=param(18,1); Le=param(19,1); n=param(20,1); thw=param(21,1); L=param(22,1); cp_gas=param(23,1); Rp=param(24,1); c1f=param(25,1); Bi=param(26,1); dPdt=param(27,1); bp1a=param(28,1); bp2a=param(29,1); bp1b=param(30,1); bp2b=param(31,1); Tf=param(32,1); qm1a=param(33,1); qm2a=param(34,1); qm1b=param(35,1); qm2b=param(36,1); qe1f=param(37,1); qe2f=param(38,1); qe10=param(39,1); qe20=param(40,1); A3=param(41,1); velo_f=param(42,1); P=param(43,1); cf=param(44,1); R=param(45,1);
const1=((1-eb)/eb)*tr; const2=1/(Peth*(1+Le)); const3=(1-(1/(1+Le))); dPdt=0; const4=(tr/c1f)*rho_particle*qe1f; const5=(tr/(1-c1f))*rho_particle*qe20; const6=((1-eb)/(eb*(1+Le)))*(1/(rho_gas*cp_gas*Tf)); %------------------------------------------------------------- dudt=zeros(7*n,1); dc1dt=zeros(7*n,1); dvelodt=zeros(7*n,1); dcp1dt=zeros(7*n,1); dcp2dt=zeros(7*n,1); dq1dt=zeros(7*n,1); dq2dt=zeros(7*n,1); dTdt=zeros(7*n,1); %--------------------------------------------------------------- c1(1:n)=u(1:n); velo(1:n)=u(n+1:2*n); cp1(1:n)=u(2*n+1:3*n); cp2(1:n)=u(3*n+1:4*n); q1(1:n)=u(4*n+1:5*n); q2(1:n)=u(5*n+1:6*n); T(1:n)=u(6*n+1:7*n); %------------------------------------------------------- %qe1 = (qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp2(i) ); %qe2 = (qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp2(i) ); %------------------------------------------------------- for i=1:n if i==1 %BOUNDARY CONDITIONS AT X=0 %--------------------------------------------------------- dc1dt(i)=(1/Pe)*((c1(i+1)-2*c1(i)+c1f)/(dx^2))... -velo(i)*((c1(i+1)-c1f)/2*dx)... -c1(i)*((velo(i+1)-velo_f)/2*dx)... +const1*k1*(c1(i)-cp1(i));
dvelodt(i)=const2*(T(i+1)-2*T(i)+Tf)/(dx^2)...
+(2/(Pe*T(i)))*(T(i+1)+Tf)/(2*dx)...
-(1/Pe)*(T(i+1)-2*T(i)+Tf)/(dx^2) ...
+velo(i)*const3*(T(i+1)-Tf)/(2*dx)...
-(T(i)/P)*dPdt...
-T(i)*(velo(i+1)-velo_f)/(2*dx)...
-A1*dq1dt(i)...
-A3*dq2dt(i)...
+A2*(thw-T(i))...
+F*(T(i)^2)*k1*tr*(c1(i)-cp1(i))...
+F*(T(i)^2)*k2*tr*((1-c1(i))-cp2(i));
dcp1dt(i)=((k1*tr)/ep)*(c1(i)-cp1(i))
-kc1*const4*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dcp2dt(i)=((k2*tr)/ep)*((1-c1(i))-cp1(i))
-kc2*const5*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dq1dt(i)=kc1*tr(i)*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dq2dt(i)=kc2*tr*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dTdt(i)=const2*(T(i+1)-2*T(i)+Tf)/(dx^2)...
-(velo(i)/(1+Le))*(T(i+1)-Tf)/(2*dx)...
-const6*kc1*tr*qe1f*DH1*(((qm1a+qm1b*(Tf*T(i)))*...
bp1a*exp(bp1b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp1(i) /...
(1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp1(i)+...
(bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i))...
-const6*kc2*tr*qe20*DH2*(((qm2a+qm2b*(Tf*T(i)))*...
bp2a*exp(bp2b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp2(i) /...
( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i) )/qe20)-q2(i))
+A2*(thw-T(i));
%---------------------------------------------------------
elseif i==2:n-1
%--------------------------------------------------------
dc1dt(i)=(1/Pe)*((c1(i+1)-2*c1(i)+c1(i-1))/(dx^2))...
-velo(i)*((c1(i+1)-c1(i-1))/2*dx)...
-c1(i)*((velo(i+1)-velo(i-1))/2*dx)...
+const1*k1*(c1(i)-cp1(i));
dvelodt(i)=const2*(T(i+1)-2*T(i)+T(i-1))/(dx^2)...
+(2/(Pe*T(i)))*(T(i+1)+T(i-1))/(2*dx)...
-(1/Pe)*(T(i+1)-2*T(i)+T(i-1))/(dx^2) ...
+velo(i)*const3*(T(i+1)-T(i-1))/(2*dx)...
-(T(i)/P)*dPdt...
-T(i)*(velo(i+1)-velo(i-1))/(2*dx)...
-A1*dq1dt(i)...
-A3*dq2dt(i)...
+A2*(Tw-T(i))...
+F*(T(i)^2)*k1*tr*(c1(i)-cp1(i))...
+F*(T(i)^2)*k2*tr*((1-c1(i))-cp2(i));
dcp1dt(i)=((k1*tr)/ep)*(c1(i)-cp1(i))
-kc1*const4*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dcp2dt(i)=((k2*tr)/ep)*((1-c1(i))-cp1(i))
-kc2*const5*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i))/qe20)-q2(i));
dq1dt(i)=kc1*tr(i)*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dq2dt(i)=kc2*tr*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dTdt(i)=const2*(T(i+1)-2*T(i)+T(i-1))/(dx^2)...
-(velo(i)/(1+Le))*(T(i+1)-T(i-1))/(2*dx)...
-const6*kc1*tr*qe1f*DH1*(((qm1a+qm1b*(Tf*T(i)))*...
bp1a*exp(bp1b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp1(i)...
/( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i))/qe1f)-q1(i))
-const6*kc2*tr*qe20*DH2*(((qm2a+qm2b*(Tf*T(i)))*...
bp2a*exp(bp2b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp2(i) /...
( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i))/qe20)-q2(i))
+A2*(thw-T(i));
elseif i==n
%BOUNDARY CONDITIONS AT X=L
%--------------------------------------------------------
dc1dt(i)=(1/Pe)*((-2*c1(i)+2*c1(i-1))/(dx^2))...
+const1*k1*(c1(i)-cp1(i));
dvelodt(i)=const2*(-2*T(i)+2*T(i-1))/(dx^2)...
-(1/Pe)*(-2*T(i)+T(i-1))/(dx^2) ...
-(T(i)/P)*dPdt...
-A1*dq1dt(i)...
-A3*dq2dt(i)...
+A2*(thw-T(i))...
+F*(T(i)^2)*k1*tr*(c1(i)-cp1(i))...
+F*(T(i)^2)*k2*tr*((1-c1(i))-cp2(i));
dcp1dt(i)=((k1*tr)/ep)*(c1(i)-cp1(i))
-kc1*const4*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dcp2dt(i)= ((k2*tr)/ep)*((1-c1(i))-cp1(i))
-kc2*const5*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dq1dt(i)=kc1*tr*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dq2dt(i)=kc2*tr*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dTdt(i)=const2*(-2*T(i)+2*T(i-1))/(dx^2)...
-const6*kc1*tr*qe1f*DH1*(((qm1a+qm1b*(Tf*T(i)))*...
bp1a*exp(bp1b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp1(i) /...
( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i) )/qe1f)-q1(i))...
-const6*kc2*tr*qe20*DH2*...
(((qm2a+qm2b*(Tf*T(i)))*...
bp2a*exp(bp2b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp2(i) /...
( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i) )/qe20)-q2(i))...
+A2*(thw-T(i));
%--------------------------------------------------------
end
dudt(i)=dc1dt(i);
dudt(n+i)=dvelodt(i);
dudt(2*n+i)=dcp1dt(i);
dudt(3*n+i)=dcp2dt(i);
dudt(4*n+i)=dq1dt(i);
dudt(5*n+i)=dq2dt(i);
dudt(6*n+i)=dTdt(i);
end %-------------------------------------------------------------
end
end
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/193759/image.png)
1 Comment
Jan
on 14 Aug 2018
Please use the "{} Code" button to format the code. Post the complete error message, such that the readers can know in which line of which function the problem occurs.
Note: 10^5 is a very expensive power operation, while 1e5 is a cheap constant.
Accepted Answer
More Answers (1)
Jan
on 14 Aug 2018
Although the code is hard to read yet, the error "Unable to meet integration tolerances" means most likely that the function to be integrated has a pole or a discontinuity at this point. Simply try this by integrating until 1.1306e-05 and plot the results. What do you observe?
0 Comments
See Also
Categories
Find more on Particle & Nuclear Physics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!