第6章_偏微分方程数值解法(附录)

1第6章附录6.1对流方程例题6.1.2计算程序!ADVECT.F!!!!advect-Programtosolvetheadvectionequation!usingthevarioushyperbolicPDEschemesprogramadvectinteger*4MAXN,MAXnplotsparameter(MAXN=500,MAXnplots=500)integer*4method,N,nStep,i,j,ip(MAXN),im(MAXN)integer*4iplot,nplots,iStepreal*8L,h,c,tau,coeff,coefflw,pi,sigma,k_wavereal*8x(MAXN),a(MAXN),a_new(MAXN),plotStepreal*8aplot(MAXN,MAXnplots+1),tplot(MAXnplots+1)!*Selectnumericalparameters(timestep,gridspacing,etc.).write(*,*)'Chooseanumericalmethod:'write(*,*)'1)FTCS,2)Lax,3)Lax-Wendroff:'read(*,*)methodwrite(*,*)'Enternumberofgridpoints:100'read(*,*)NL=1.!Systemsizeh=L/N!Gridspacingc=1!Wavespeedwrite(*,*)'Timeforwavetomoveonegridspacingis',h/cwrite(*,*)'Entertimestep:0.001'read(*,*)taucoeff=-c*tau/(2.*h)!Coefficientusedbyallschemescoefflw=2*coeff*coeff!CoefficientusedbyL-Wschemewrite(*,*)'Wavecirclessystemin',L/(c*tau),'steps'write(*,*)'Enternumberofsteps:300'read(*,*)nStep!*Setinitialandboundaryconditions.pi=3.141592654sigma=0.1!WidthoftheGaussianpulsek_wave=pi/sigma!Wavenumberofthecosinedoi=1,Nx(i)=(i-0.5)*h-L/2!Coordinatesofgridpoints!InitialconditionisaGaussian-cosinepulse2a(i)=cos(k_wave*x(i))*exp(-x(i)**2/(2*sigma**2))enddo!Useperiodicboundaryconditionsdoi=2,(N-1)ip(i)=i+1!ip(i)=i+1withperiodicb.c.im(i)=i-1!im(i)=i-1withperiodicb.c.enddoip(1)=2ip(N)=1im(1)=Nim(N)=N-1!*Initializeplottingvariables.iplot=1!Plotcounternplots=50!DesirednumberofplotsplotStep=float(nStep)/nplotstplot(1)=0!Recordtheinitialtime(t=0)doi=1,Naplot(i,1)=a(i)!Recordtheinitialstateenddo!*Loopoverdesirednumberofsteps.doiStep=1,nStep!*ComputenewvaluesofwaveamplitudeusingFTCS,!LaxorLax-Wendroffmethod.if(method.eq.1)then!!!FTCSmethod!!!doi=1,Na_new(i)=a(i)+coeff*(a(ip(i))-a(im(i)))enddoelseif(method.eq.2)then!!!Laxmethod!!!doi=1,Na_new(i)=0.5*(a(ip(i))+a(im(i)))+&coeff*(a(ip(i))-a(im(i)))enddoelse!!!Lax-Wendroffmethod!!!doi=1,Na_new(i)=a(i)+coeff*(a(ip(i))-a(im(i)))+&coefflw*(a(ip(i))+a(im(i))-2*a(i))enddoendifdoi=1,N3a(i)=a_new(i)!Resetwithnewamplitudevaluesenddo!*Periodicallyrecorda(t)forplotting.if((iStep-int(iStep/plotStep)*plotStep).lt.1)theniplot=iplot+1tplot(iplot)=tau*iStepdoi=1,Naplot(i,iplot)=a(i)!Recorda(i)forplotingenddowrite(*,*)iStep,'outof',nStep,'stepscompleted'endifenddonplots=iplot!Actualnumberofplotsrecorded!*Printouttheplottingvariables:x,a,tplot,aplotopen(11,file='x.txt',status='unknown')open(12,file='a.txt',status='unknown')open(13,file='tplot.txt',status='unknown')open(14,file='aplot.txt',status='unknown')doi=1,Nwrite(11,*)x(i)write(12,*)a(i)doj=1,(nplots-1)write(14,1001)aplot(i,j)enddowrite(14,*)aplot(i,nplots)enddodoi=1,nplotswrite(13,*)tplot(i)enddo1001format(e12.6,',',$)!The$suppressesthecarriagereturnstopend!*****ToplotinMATLAB;usethescriptbelow********************!loadx.txt;loada.txt;loadtplot.txt;loadaplot.txt;!%*Plottheinitialandfinalstates.!figure(1);clf;%Clearfigure1windowandbringforward!plot(x,aplot(:,1),'-',x,a,'--');!legend('Initial','Final');!xlabel('x');ylabel('a(x,t)');4!pause(1);%Pause1secondbetweenplots!%*Plotthewaveamplitudeversuspositionandtime!figure(2);clf;%Clearfigure2windowandbringforward!mesh(tplot,x,aplot);!ylabel('Position');xlabel('Time');zlabel('Amplitude');!view([-7050]);%Betterviewfromthisangle!******************************************************************================================================6.2抛物形方程例题6.2.2计算程序!!!!dftcs.f!!!!!!!!!!!dftcs-Programtosolvethediffusionequation!usingtheForwardTimeCenteredSpace(FTCS)scheme.programdftcsinteger*4MAXN,MAXnplotsparameter(MAXN=300,MAXnplots=500)integer*4N,i,j,iplot,nStep,plot_step,nplots,iStepreal*8tau,L,h,kappa,coeff,tt(MAXN),tt_new(MAXN)real*8xplot(MAXN),tplot(MAXnplots),ttplot(MAXN,MAXnplots)!*Initializeparameters(timestep,gridspacing,etc.).write(*,*)'Entertimestep:0.0001'read(*,*)tauwrite(*,*)'Enterthenumberofgridpoints:50'read(*,*)NL=1.!Thesystemextendsfromx=-L/2tox=L/2h=L/(N-1)!Gridsizekappa=1.!Diffusioncoefficientcoeff=kappa*tau/h**2if(coeff.lt.0.5)thenwrite(*,*)'Solutionisexpectedtobestable'elsewrite(*,*)'WARNING:Solutionisexpectedtobeunstable'endif!*Setinitialandboundaryconditions.doi=1,Ntt(i)=0.0!Initializetemperaturetozeroatallpointstt_new(i)=0.0enddott(N/2)=1/h!Initialcond.isdeltafunctionincenter5!!Theboundaryconditionsarett(1)=tt(N)=0!!Endpointsareunchangedduringiteration!*Setuploopandplotvariables.iplot=1!CounterusedtocountplotsnStep=300!Maximumnumberofiterationsplot_step=6!Numberoftimestepsbetweenplotsnplots=nStep/plot_step+1!Numberofsnapshots(plots)doi=1,Nxplot(i)=(i-1)*h-L/2!Recordthexscaleforplotsenddo!*Loopoverthedesirednumberoftimesteps.doiStep=1,nStep!*ComputenewtemperatureusingFTCSscheme.doi=2,(N-1)tt_new(i)=tt(i)+coeff*(tt(i+1)+tt(i-1)-2*tt(i))enddodoi=2,(N-1)tt(i)=tt_new(i)!Resettemperaturetonewvaluesenddo!*Periodicallyrecordtemperatureforplotting.if(mod(iStep,plot_step).lt.1)then!Everyplot_stepstepsdoi=1,N!recordtt(i)forplottingttplot(i,iplot)=tt(i)enddotplot(iplot)=iStep*tau!Recordtimeforplotsiplot=iplot+1endifenddonplots=iplot-1!Numberofplotsactuallyrecorded!*Printouttheplottin