ODE Solvers fail to converge as the tolerances are decreased.

Question

Zachary Epstein 2017-8-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/354509-ode-solvers-fail-to-converge-as-the-tolerances-are-decreased

评论： Zachary Epstein 2017-8-30

I'm numerically integrating a coupled ODE with time-dependent coefficients. When the ODE splits the time window into a good number of steps (say, MaxStep=0.5), everything begins to look okay - but as the time steps get smaller (say, MaxStep=0.01) the solution starts looking bad instead of converging. My theory is that the function is driven by small numbers which border Matlab's digit limit, such that for tiny time steps the function takes longer than it should to turn on.. but I don't know how to get around such a problem and perform an accurate numerical integration.

I will attach the code text below as well as my integration script, shortIntegralScript.m in case you wish to look at the code and/or run it. You will see plots which compare the numerical result to the true result. I should note that in the absence of the MaxStep option, only ode23 yields results which look reasonable, and they get worse as AbsTol increases.

Thank you!

%% Just defining the vector 'gamma'.

LowLine = 0.2482; DeltaLambda = 0.0004; c=3.*10^10; w0=2.*pi*c/((LowLine+DeltaLambda/2)*10^-4); w0=w0*10^-12; T = 150*2*pi; nt = 2^13; dt = T/nt; % timestep (dt) t = ((1:nt)'-(nt+1)/2)*dt; % time vector w = wspace(T,nt); % angular frequency vector vs = fftshift(w/(2*pi)); z = 10200.0/2; nz = 2^2; dz = z/nz; scal = 1; Aeff=5^2; Pressure = 760/5; %Torr n0 = 1+(.00084)*(Pressure/760); A0 = 0.77*sqrt(225/Aeff)*100762./(scal); Pavg=(A0/57861.7)^2*n0*Aeff/225; Norm=A0/sqrt(5); XDL=1; rcoh= 0.664*(sqrt(Aeff)/XDL); % Use this for ArrayuG LinewidthX = 2*pi*(0.664/rcoh); gmX = LinewidthX/2.35; LinewidthY = 2*pi*(0.664/rcoh); gmY = LinewidthY/2.35; LinewidthXLor = 2.*pi*((1/pi)/rcoh); gmXLor=LinewidthXLor/2; LinewidthYLor = 2.*pi*((1/pi)/rcoh); gmYLor=LinewidthYLor/2; rx = sqrt(Aeff)/2; %cm ry = sqrt(Aeff)/2; %cm Room=1.5; nx=2*ceil((XDL/2)*Room); ny=2*ceil((XDL/2)*Room); %number of points on each side of the axis X=rx*Room; Y=ry*Room; %in cm, on each side of the axis wn = w0+fftshift((transpose([-nt/2+1/2:nt/2-1/2])*(2*pi/T))); nn=[-nt/2+1/2:nt/2-1/2]'; x = ((1:(2*nx+1))'-(nx+1))*(X/nx); kxn = 0+(transpose([-nx:nx])*(2*pi/(2*X))); nnx=[-nx:nx]'; y = ( (1:(2*ny+1))'-(ny+1))*(Y/ny); kyn = 0+(transpose([-ny:ny])*(2*pi/(2*Y))); nny=[-ny:ny]'; Normalize = 500; Linewidth = 2*pi*1.00; gm = Linewidth/2.35; gmLor = Linewidth/2;%i.e. \Delta\omega in THz tcoh0G=2*pi*0.664/Linewidth; tcoh0L=2*pi*(1/pi)/Linewidth; NOscGAll=zeros(nt,2*nx+1,2*ny+1); N0iseGAll=zeros(nt,2*nx+1,2*ny+1); for ikx=1:2*nx+1 for iky=1:2*ny+1 for iw=1:nt NOscGAll(iw,ikx,iky) = (1/sqrt(2*pi*gm^2))*exp(-(1/2)*((wn(iw)-w0)./gm).^2)*(1/sqrt(2*pi*gmX^2))*exp(-(1/2)*((kxn(ikx)-0)/gmX).^2)*(1/sqrt(2*pi*gmY^2))*exp(-(1/2)*((kyn(iky)-0)/gmY).^2); end end end UTildeGAll = normrnd(0,NOscGAll).*exp(1i*2*pi*rand(nt,2*nx+1,2*ny+1)); N0iseGAll=fftn(UTildeGAll); NOscGAll=[];UTildeGAll=[]; N0iseGAll=N0iseGAll./sqrt(sum(sum(sum((abs(N0iseGAll).^2))))/(nt*(2*nx+1)*(2*ny+1))); FilterDepth=13; rt = (1/2)*T/1.5; FilterX=((1/2)*(tanh((1/rx)*FilterDepth*(x+rx))+tanh(-(1/rx)*FilterDepth*(x-rx)))); FilterY=((1/2)*(tanh((1/ry)*FilterDepth*(y+ry))+tanh(-(1/ry)*FilterDepth*(y-ry)))); FilterT=((1/2)*(tanh((1/rt)*FilterDepth*(t+rt))+tanh(-(1/rt)*FilterDepth*(t-rt)))); Filter=zeros(nt,2*nx+1,2*ny+1); for ix = 1:(2*nx+1) for iy = 1:(2*ny+1) for iz = 1:nt Filter(iz,ix,iy)=FilterX(ix)*FilterY(iy)*FilterT(iz); end end end FilterT=[];FilterX=[];FilterY=[]; Arrayu0=A0*(N0iseGAll).*Filter; Arrayu0=A0*Filter; FilterDepth=13/2; om21=2*(2*pi*3*10^5/(248.4))*0.849; om31=2*2*pi*3*10^5/(249.6); om32=om31-om21; f1=1; mu12=14.53*10^-30*100*f1^(1/4); mu23=mu12; px=(nx+1); py=(nx+1); u0=Arrayu0(:,px,py); Dv=angle(u0); t=((1:nt)'-(nt+1)/2)*dt; Dvsp=Dv-Dv(1); Numb=0; for it=2:nt if Dv(it)-Dv(it-1)>5 Numb=Numb+1; tau=[tau t(it)]; end if Dv(it)-Dv(it-1)<-5 Numb=Numb-1; tau=[tau t(it)]; end Dvsp(it)=Dvsp(it)-Numb*2*pi; end Phi=Dvsp; %(really should be + w0 t) dPhidt=w0+(Dvsp-circshift(Dvsp,1))/dt; dPhidt(1)=0; Eabs=abs(u0)/2; dPhidt(1)=dPhidt(2); a=1;b=1; hbar=(6.63*10^-34)/(2*pi); Om12 = 10^-12*mu12*Eabs/hbar; %(positive phi argument). Eabs in V/cm, so Om12 is in ps^-1 now Om23 = 10^-12*(mu23)*Eabs/hbar; %(positive phi argument) Om12dE = 10^-12*mu12/hbar; %(positive phi argument). Eabs in V/cm, so Om12 is in ps^-1 now Om23dE = 10^-12*(mu23)/hbar; %(positive phi argument) gam1=2*Om12.*Om23./((om21-dPhidt)); gam2=0*t.^0; dw31=abs(Om12).^2.*(1./(dPhidt-om32)+1./(3*dPhidt-om32))+abs(Om23).^2.*(1./(dPhidt-om21)+1./(3*dPhidt-om21)); gam3=-(om31-2*dPhidt+dw31); r1p=-sign(gam3).*gam1./sqrt(gam1.^2+gam3.^2); r3p=-sign(gam3).*gam3./sqrt(gam1.^2+gam3.^2); r2p=-((1/dt)./gam3).*(r1p-circshift(r1p,1)); r3pp=-sqrt(1-r1p.^2); %% Defining the time-dependent coefficients of the coupled ODE then Integrating N11=0*gam1;N1d=N11;N2d=N11;N3d=N11; N12=-gam3; N13=gam2; N21=gam3; N22=N11; N23=-gam1; N31=-gam2; N32=gam1; N33=N11; N11(1)=N11(2);N12(1)=N12(2);N13(1)=N13(2); N21(1)=N21(2);

r=shortIntegralScript(t,N11,N12,N13,N21,N22,N23,N31,N32,N33,0,0,-1);

r3=r(:,3); r1=r(:,1); r2=r(:,2);

r1=interp1((-length(r1)/2+1/2 : length(r1)/2-1/2)*T/length(r1),r1,t); r2=interp1((-length(r2)/2+1/2 : length(r2)/2-1/2)*T/length(r2),r2,t); r3=interp1((-length(r3)/2+1/2 : length(r3)/2-1/2)*T/length(r3),r3,t);

r1=real(r1);r2=real(r2);r3=real(r3);

figure() hold on plot(t,r1) plot(t,r1p) title(['r2: Should match']) hold off

figure() hold on plot(t,r2) plot(t,r2p) title(['r2: Should match']) hold off

figure() hold on plot(t,r3) plot(t,r3p) title(['r3: Should match']) hold off

function r = shortIntegralScript(tvar,N11,N12,N13,N21,N22,N23,N31,N32,N33,r10,r20,r30); tic nt=length(N11); T=(nt/(nt-1))*(max(tvar)-min(tvar)); dt=T/nt;

nt=length(tvar); T=(nt/(nt-1))*(max(tvar)-min(tvar)); t0=min(tvar); Methd=2; %A fast interpolation method function drdt = myode(tt,r,T,t0,nt,tvar,N11,N12,N13,N21,N22,N23,N31,N32,N33)

if Methd==0

aN11 = interp1(tvar,N11,tt,'spline'); % Interpolate the data set (ft,f) at time t aN12 = interp1(tvar,N12,tt,'spline'); % Interpolate the data set (ft,f) at time t aN13 = interp1(tvar,N13,tt,'spline'); % Interpolate the data set (ft,f) at time t aN21 = interp1(tvar,N21,tt,'spline'); % Interpolate the data set (ft,f) at time t aN22 = interp1(tvar,N22,tt,'spline'); % Interpolate the data set (ft,f) at time t aN23 = interp1(tvar,N23,tt,'spline'); % Interpolate the data set (ft,f) at time t aN31 = interp1(tvar,N31,tt,'spline'); % Interpolate the data set (ft,f) at time t aN32 = interp1(tvar,N32,tt,'spline'); % Interpolate the data set (ft,f) at time t aN33 = interp1(tvar,N33,tt,'spline'); % Interpolate the data set (ft,f) at time t

elseif Methd==1 aN11 = interp1(tvar,N11,tt); % Interpolate the data set (ft,f) at time t aN12 = interp1(tvar,N12,tt); % Interpolate the data set (ft,f) at time t aN13 = interp1(tvar,N13,tt); % Interpolate the data set (ft,f) at time t aN21 = interp1(tvar,N21,tt); % Interpolate the data set (ft,f) at time t aN22 = interp1(tvar,N22,tt); % Interpolate the data set (ft,f) at time t aN23 = interp1(tvar,N23,tt); % Interpolate the data set (ft,f) at time t aN31 = interp1(tvar,N31,tt); % Interpolate the data set (ft,f) at time t aN32 = interp1(tvar,N32,tt); % Interpolate the data set (ft,f) at time t aN33 = interp1(tvar,N33,tt); % Interpolate the data set (ft,f) at time t

elseif Methd==2 ValT=ceil(nt*(tt-t0)/T); Per=ceil(nt*(tt-t0)/T)-nt*(tt-t0)/T; if ValT==0 ValT=1; aN11 = N11(ValT); aN12 = N12(ValT); aN13 = N13(ValT); aN21 = N21(ValT); aN22 = N22(ValT); aN23 = N23(ValT); aN31 = N31(ValT); aN32 = N32(ValT); aN33 = N33(ValT); elseif ValT==nt aN11 = N11(ValT); aN12 = N12(ValT); aN13 = N13(ValT); aN21 = N21(ValT); aN22 = N22(ValT); aN23 = N23(ValT); aN31 = N31(ValT); aN32 = N32(ValT); aN33 = N33(ValT); else aN11 = Per*N11(ValT)+(1-Per)*N11(ValT+1); aN12 = Per*N12(ValT)+(1-Per)*N12(ValT+1); aN13 = Per*N13(ValT)+(1-Per)*N13(ValT+1); aN21 = Per*N21(ValT)+(1-Per)*N21(ValT+1); aN22 = Per*N22(ValT)+(1-Per)*N22(ValT+1); aN23 = Per*N23(ValT)+(1-Per)*N23(ValT+1); aN31 = Per*N31(ValT)+(1-Per)*N31(ValT+1); aN32 = Per*N32(ValT)+(1-Per)*N32(ValT+1); aN33 = Per*N33(ValT)+(1-Per)*N33(ValT+1); end

elseif Methd==3 ValT=ceil(nt*(tt-t0)/T); Per=ceil(nt*(tt-t0)/T)-nt*(tt-t0)/T; if ValT==0 ValT=1; end aN11 = N11(ValT); aN12 = N12(ValT); aN13 = N13(ValT); aN21 = N21(ValT); aN22 = N22(ValT); aN23 = N23(ValT); aN31 = N31(ValT); aN32 = N32(ValT); aN33 = N33(ValT); end

drdt = zeros(3,1); drdt(1) = aN11.*r(1)+aN12.*r(2)+aN13.*r(3); drdt(2) = aN21.*r(1)+aN22.*r(2)+aN23.*r(3); drdt(3) = aN31.*r(1)+aN32.*r(2)+aN33.*r(3);

end

tspan = [min(tvar) max(tvar)]; ic = [0 0 -1]; opts = odeset('RelTol', 1e-8,'AbsTol',1e-10,'MaxStep',(.01)) %Here is where things get crazy. [tt,r] = ode23(@(tt,r) myode(tt,r,T,t0,nt,tvar,N11,N12,N13,N21,N22,N23,N31,N32,N33), tspan, ic,opts);

end

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

David Goodmanson 2017-8-30

Hi Zachary,

You need to modify the question by highlighting the code and using the {}code button. It's long code even then, but if you don't do that, people are faced with the task of unentangling what you have.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Jan 2017-8-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/354509-ode-solvers-fail-to-converge-as-the-tolerances-are-decreased#answer_279659

编辑：Jan 2017-8-30

Your code is hard to read, but is looks like some of the parameters are interpolated linearly. Then the function be be integrated is not smooth and therefore it cannot be integrated by Matlab's ODE23 reliably. ceil and if might introduce additional discontinuities. See http://www.mathworks.com/matlabcentral/answers/59582#answer_72047. Driving an integrator outside its specifications is a really bad idea. I'm not surprised that the convergence does not work. But it would be even worse if it works: Then Rounding or discretization errors can dominate the "result" and then the integration is an inefficient random number generator.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Zachary Epstein 2017-8-30

在 MATLAB Online 中打开

Thank you very much for the response!

Your answer was very insightful, although it turns out that the stiffness is the source of the problem rather than the interpolation - Below I have simplified the code and use exact coefficients (no interpolation). The resulting plot of the numerical-vs-exact solution produces the same issue.

Although the code is for laser physics research, the simplified code can be thought of as:

A vector r rotates about vector gamma at 10 times per second. The vector gamma slowly changes direction on the order of 100 seconds. I want to integrate the system numerically to determine the path of r over 900 seconds.

I have looked at all the main Matlab ODE Solvers, and nothing seems to be able to resolve the stiffness. The best I can do is ode23, and even that requires me to set a MaxStep of 0.005-0.5 - and I'm concerned that that counts as 'driving an integrator outside its specifications'.

Thank you very much.

- Zach

Here I just define the rotation vector 'gamma'.

w0=7588.388; T = 900.; nt = 2^13; dt = T/nt; t = ((1:nt)'-(nt+1)/2)*dt;   A0 = 232760; FilterDepth=13; rt = T/3; om21=12885;
  om31=15104;
  om32=22187;
  mu12=14.53*10^-40;mu23=mu12;
  Eabs=abs(A0*((1/2)*(tanh((1/rt)*FilterDepth*(t+rt))+tanh(-(1/rt)*FilterDepth*(t-rt)))))/2;
  hbar=(6.63*10^-34)/(2*pi);
  Om12 = mu12*Eabs/hbar; %(positive phi argument). Eabs in V/cm, so Om12 is in ps^-1 now
  Om23 = (mu23)*Eabs/hbar; %(positive phi argument)
  dw31=abs(Om12).^2.*(1./(w0-om32)+1./(3*w0-om32))+abs(Om23).^2.*(1./(w0-om21)+1./(3*w0-om21));
gam1=2*Om12.*Om23./((om21-w0));
gam3=-(om31-2*w0+dw31);}
% Here I run the numerical integration
r=ZachshortIntegralScript(t,w0,T);
% Here I compare the numerical result to the expected solution rp
r3=r(:,3); r1=r(:,1); r2=r(:,2);
r1=interp1((-length(r1)/2+1/2 : length(r1)/2-1/2)*T/length(r1),r1,t);
r2=interp1((-length(r2)/2+1/2 : length(r2)/2-1/2)*T/length(r2),r2,t);
r3=interp1((-length(r3)/2+1/2 : length(r3)/2-1/2)*T/length(r3),r3,t);
r1p=-sign(gam3).*gam1./sqrt(gam1.^2+gam3.^2);
r3p=-sign(gam3).*gam3./sqrt(gam1.^2+gam3.^2);
r2p=-((1/dt)./gam3).*(r1p-circshift(r1p,1));
figure()
hold on
plot(t,r1)
plot(t,r1p)
title(['r2: Should match'])
hold off
figure()
hold on
plot(t,r2)
plot(t,r2p)
title(['r2: Should match'])
hold off
figure()
hold on
plot(t,r3)
plot(t,r3p)
title(['r3: Should match'])
hold off

To conclude, here is my numerical integration script

function r = ZachshortIntegralScript(tvar,w0,T);
nt=length(tvar); t0=min(tvar);
function drdt = myode(tt,r,t0,nt,tvar)
% Just defining gam1 and gam3 again.
A0 = 232760; FilterDepth=13; rt = T/3; om21=12885;
om31=15104; om32=22187;
mu12=14.53*10^-40;mu23=mu12;
Eabs=abs(A0*((1/2)*(tanh((1/rt)*FilterDepth*(tt+rt))+tanh(-(1/rt)*FilterDepth*(tt-rt)))))/2;
hbar=(6.63*10^-34)/(2*pi);
Om12 = mu12*Eabs/hbar; Om23 = (mu23)*Eabs/hbar; dw31=abs(Om12).^2.*(1./(w0-om32)+1./(3*w0-om32))+abs(Om23).^2.*(1./(w0-om21)+1./(3*w0-om21)); 
gam1=2*Om12.*Om23./((om21-w0));
gam3=-(om31-2*w0+dw31);
drdt = zeros(3,1);
drdt(1) = -gam3.*r(2);
drdt(2) = gam3.*r(1)+-gam1.*r(3);
drdt(3) = gam1.*r(2);
end
tspan = [min(tvar) max(tvar)];
ic = [0 0 -1];
opts = odeset('MaxStep',(.001))
[tt r] = ode45(@(tt,r) myode(tt,r,t0,nt,tvar), tspan, ic);
end