ODE Solvers fail to converge as the tolerances are decreased.
5 次查看(过去 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 个评论
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.
回答(1 个)
Jan
2017-8-30
编辑: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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!