偏微分方程数值解法答案
偏微分方程数值解期末试题及标准答案
偏微分方程数值解试题(０6B ）参考答案与评分标准信息与计算科学专业一(10分）、设矩阵A 对称，定义)(),(),(21)(n R x x b x Ax x J ∈-=，)()(0x x J λλϕ+=.若0)0('=ϕ，则称称0x 是)(x J 的驻点(或稳定点).矩阵A 对称(不必正定)，求证0x 是)(x J 的驻点的充要条件是:0x 是方程组 b Ax =的解 解: 设n R x ∈0是)(x J 的驻点，对于任意的n R x ∈,令),(2),()()()(2000x Ax x b Ax x J x x J λλλλϕ+-+=+=, （3分)0)0('=ϕ,即对于任意的n R x ∈,0),(0=-x b Ax ,特别取b Ax x -=0，则有0||||),(2000=-=--b Ax b Ax b Ax ，得到b Ax =0． (3分)反之,若n R x ∈0满足b Ax =0,则对于任意的x ,)(),(21)0()1()(00x J x Ax x x J >+==+ϕϕ,因此0x 是)(x J 的最小值点. (４分）评分标准：)(λϕ的展开式3分, 每问3分，推理逻辑性１分二（10分）、 对于两点边值问题:⎪⎩⎪⎨⎧==∈=+-=0)(,0)(),()('b u a u b a x f qu dx du p dx d Lu 其中]),([,0]),,([,0)(min )(]),,([0min ],[1b a H f q b a C q p x p x p b a C p b a x ∈≥∈>=≥∈∈建立与上述两点边值问题等价的变分问题的两种形式：求泛函极小的Ritz 形式和Galerkin 形式的变分方程。
解: 设}0)(),,(|{11=∈=a u b a H u u H E为求解函数空间,检验函数空间．取),(1b a H v E ∈,乘方程两端,积分应用分部积分得到 (３分))().(),(v f fvdx dx quv dxdv dx du p v u a b a ba ==+=⎰⎰,),(1b a H v E ∈∀ 即变分问题的Galerkin 形式. （3分)令⎰-+=-=b a dx fu qu dxdu p u f u u a u J ])([21),(),(21)(22，则变分问题的Ritz 形式为求),(1*b a H u E ∈,使)(m in )(1*u J u J EH u ∈= （4分) 评分标准:空间描述与积分步骤３分，变分方程3分，极小函数及其变分问题４分,三(20分)、对于边值问题⎪⎩⎪⎨⎧-====⨯=∈=∂∂+∂∂====x u u u u G y x y u x u y y x x 1||,0|,1|)1,0()1,0(),(,010102222 (１)建立该边值问题的五点差分格式(五点棱形格式又称正五点格式),推导截断误差的阶。
第6章_偏微分方程数值解法(附录)
第6章附录6.1 对流方程例题6.1.2计算程序!ADVECT.F!!!! advect - Program to solve the advection equation! using the various hyperbolic PDE schemesprogram advectinteger*4 MAXN, MAXnplotsparameter( MAXN = 500, MAXnplots = 500 )integer*4 method, N, nStep, i, j, ip(MAXN), im(MAXN)integer*4 iplot, nplots, iStepreal*8 L, h, c, tau, coeff, coefflw, pi, sigma, k_wavereal*8 x(MAXN), a(MAXN), a_new(MAXN), plotStepreal*8 aplot(MAXN,MAXnplots+1), tplot(MAXnplots+1)!* Select numerical parameters (time step, grid spacing, etc.).write(*,*) 'Choose a numerical method: 'write(*,*) ' 1) FTCS, 2) Lax, 3) Lax-Wendroff : 'read(*,*) methodwrite(*,*) 'Enter number of grid points: 100'read(*,*) NL = 1. ! System sizeh = L/N ! Grid spacingc = 1 ! Wave speedwrite(*,*) 'Time for wave to move one grid spacing is ', h/cwrite(*,*) 'Enter time step: 0.001'read(*,*) taucoeff = -c*tau/(2.*h) ! Coefficient used by all schemescoefflw = 2*coeff*coeff ! Coefficient used by L-W schemewrite(*,*) 'Wave circles system in ', L/(c*tau), ' steps'write(*,*) 'Enter number of steps: 300'read(*,*) nStep!* Set initial and boundary conditions.pi = 3.141592654sigma = 0.1 ! Width of the Gaussian pulsek_wave = pi/sigma ! Wave number of the cosinedo i=1,Nx(i) = (i-0.5)*h - L/2 ! Coordinates of grid pointsa(i) = cos(k_wave*x(i)) * exp(-x(i)**2/(2*sigma**2))enddo! Use periodic boundary conditionsdo i=2,(N-1)ip(i) = i+1 ! ip(i) = i+1 with periodic b.c.im(i) = i-1 ! im(i) = i-1 with periodic b.c.enddoip(1) = 2ip(N) = 1im(1) = Nim(N) = N-1!* Initialize plotting variables.iplot = 1 ! Plot counternplots = 50 ! Desired number of plotsplotStep = float(nStep)/nplotstplot(1) = 0 ! Record the initial time (t=0)do i=1,Naplot(i,1) = a(i) ! Record the initial stateenddo!* Loop over desired number of steps.do iStep=1,nStep!* Compute new values of wave amplitude using FTCS,! Lax or Lax-Wendroff method.if( method .eq. 1 ) then !!! FTCS method !!!do i=1,Na_new(i) = a(i) + coeff*( a(ip(i))-a(im(i)) )enddoelse if( method .eq. 2 ) then !!! Lax method !!!do i=1,Na_new(i) = 0.5*( a(ip(i))+a(im(i)) ) +& coeff*( a(ip(i))-a(im(i)) )enddoelse !!! Lax-Wendroff method !!!do i=1,Na_new(i) = a(i) + coeff*( a(ip(i))-a(im(i)) ) +& coefflw*( a(ip(i))+a(im(i))-2*a(i) ) enddoendifa(i) = a_new(i) ! Reset with new amplitude valuesenddo!* Periodically record a(t) for plotting.if( (iStep-int(iStep/plotStep)*plotStep) .lt. 1 ) theniplot = iplot+1tplot(iplot) = tau*iStepdo i=1,Naplot(i,iplot) = a(i) ! Record a(i) for plotingenddowrite(*,*) iStep, ' out of ', nStep, ' steps completed'endifenddonplots = iplot ! Actual number of plots recorded!* Print out the plotting variables: 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')do i=1,Nwrite(11,*) x(i)write(12,*) a(i)do j=1,(nplots-1)write(14,1001) aplot(i,j)enddowrite(14,*) aplot(i,nplots)enddodo i=1,nplotswrite(13,*) tplot(i)enddo1001 format(e12.6,', ',$) ! The $ suppresses the carriage return stopend!***** To plot in MATLAB; use the script below ******************** !load x.txt; load a.txt; load tplot.txt; load aplot.txt;!%* Plot the initial and final states.!figure(1); clf; % Clear figure 1 window and bring forward!plot(x,aplot(:,1),'-',x,a,'--');!legend('Initial','Final');!pause(1); % Pause 1 second between plots!%* Plot the wave amplitude versus position and time!figure(2); clf; % Clear figure 2 window and bring forward!mesh(tplot,x,aplot);!ylabel('Position'); xlabel('Time'); zlabel('Amplitude');!view([-70 50]); % Better view from this angle!****************************************************************** ================================================6.2 抛物形方程例题6.2.2计算程序!!!!dftcs.f!!!!!!!!!!! dftcs - Program to solve the diffusion equation! using the Forward Time Centered Space (FTCS) scheme.program dftcsinteger*4 MAXN, MAXnplotsparameter( MAXN = 300, MAXnplots = 500 )integer*4 N, i, j, iplot, nStep, plot_step, nplots, iStepreal*8 tau, L, h, kappa, coeff, tt(MAXN), tt_new(MAXN)real*8 xplot(MAXN), tplot(MAXnplots), ttplot(MAXN,MAXnplots)!* Initialize parameters (time step, grid spacing, etc.).write(*,*) 'Enter time step: 0.0001'read(*,*) tauwrite(*,*) 'Enter the number of grid points: 50 'read(*,*) NL = 1. ! The system extends from x=-L/2 to x=L/2h = L/(N-1) ! Grid sizekappa = 1. ! Diffusion coefficientcoeff = kappa*tau/h**2if( coeff .lt. 0.5 ) thenwrite(*,*) 'Solution is expected to be stable'elsewrite(*,*) 'WARNING: Solution is expected to be unstable'endif!* Set initial and boundary conditions.do i=1,Ntt(i) = 0.0 ! Initialize temperature to zero at all pointstt_new(i) = 0.0enddo!! The boundary conditions are tt(1) = tt(N) = 0!! End points are unchanged during iteration!* Set up loop and plot variables.iplot = 1 ! Counter used to count plotsnStep = 300 ! Maximum number of iterationsplot_step = 6 ! Number of time steps between plots nplots = nStep/plot_step + 1 ! Number of snapshots (plots)do i=1,Nxplot(i) = (i-1)*h - L/2 ! Record the x scale for plotsenddo!* Loop over the desired number of time steps.do iStep=1,nStep!* Compute new temperature using FTCS scheme.do i=2,(N-1)tt_new(i) = tt(i) + coeff*(tt(i+1) + tt(i-1) - 2*tt(i))enddodo i=2,(N-1)tt(i) = tt_new(i) ! Reset temperature to new valuesenddo!* Periodically record temperature for plotting.if( mod(iStep,plot_step) .lt. 1 ) then ! Every plot_step stepsdo i=1,N ! record tt(i) for plottingttplot(i,iplot) = tt(i)enddotplot(iplot) = iStep*tau ! Record time for plotsiplot = iplot+1endifenddonplots = iplot-1 ! Number of plots actually recorded!* Print out the plotting variables: tplot, xplot, ttplotopen(11,file='tplot.txt',status='unknown')open(12,file='xplot.txt',status='unknown')open(13,file='ttplot.txt',status='unknown')do i=1,nplotswrite(11,*) tplot(i)enddowrite(12,*) xplot(i)do j=1,(nplots-1)write(13,1001) ttplot(i,j)enddowrite(13,*) ttplot(i,nplots)enddo1001 format(e12.6,', ',$) ! The $ suppresses the carriage returnstopend!***** To plot in MATLAB; use the script below ********************!load tplot.txt; load xplot.txt; load ttplot.txt;!%* Plot temperature versus x and t as wire-mesh and contour plots.!figure(1); clf;!mesh(tplot,xplot,ttplot); % Wire-mesh surface plot!xlabel('Time'); ylabel('x'); zlabel('T(x,t)');!title('Diffusion of a delta spike');!pause(1);!figure(2); clf;!contourLevels = 0:0.5:10; contourLabels = 0:5;!cs = contour(tplot,xplot,ttplot,contourLevels); % Contour plot!clabel(cs,contourLabels); % Add labels to selected contour levels!xlabel('Time'); ylabel('x');!title('Temperature contour plot');!******************************************************************! program diffu_explicit.for! explicit methods for finite diffrencedimension T(51,51)real dx,dt,lamda,kp,Tmax,Tminopen(2,file='diffu_1.dat')dx=0.5dt=0.1kp=0.835lamda=kp*dt/dx/dxTmax=100.Tmin=0.imax=21write(*,*) 'lamda=',lamdado i=1,50do j=1,50T(i,j)=0.0do j=1,40T(1,j)=TmaxT(imax,j)=Tminend dodo j=1,40do i=2,20T(i,j+1)=T(i,j)+lamda*(T(i+1,j)-2.*T(i,j)+T(i-1,j)) end dowrite(2,222) (T(i,j),i=1,21)end do222 format(1x,21(2x,f18.10))end% diffu_explicit.mclc; clear all;nl = 20.; nt = 40;dx = 0.5; dt = 0.1; kp = 0.835;lamda = kp*dt/dx/dx;Tl = 0.; Tr = 100.;T = zeros(nl,nt);T(1,:) = Tl; T(nl,:)=Tr;for j=1:nt-1for i =2:nl-1T(i,j+1)=T(i,j)+lamda*(T(i+1,j)-2.*T(i,j)+T(i-1,j));endendsurf(T');xlabel('x');ylabel('time');zlabel('Temperature');title('explicit scheme');! diffu_implicit.for! implicit methods for finite diffrenceexternal tridparameter (imax=21,kmax=50)dimension t(imax),a(imax),b(imax),c(imax),r(imax),u(imax) real dx,dt,lamda,kp,tmax,tminopen(2,file='diffu_2.dat')dx=0.5lamda=kp*dt/dx/dxtmax=100.tmin=0.write(*,*) 'lamda=',lamdado i=2,imax-1! t(i)=tmax-(tmax-tmin)/(imax-1.)*(i-1.) t(i)=0.0end dot(1)=tmaxt(imax)=tmin! write(*,*) a(i),b(i),c(i)b(1)=1.c(1)=0.r(1)=t(1)a(imax)=0.b(imax)=1.r(imax)=t(imax)do i=2,imax-1a(i)=-lamdab(i)=1.+2.*lamdac(i)=-lamdar(i)=t(i)end dok=01 k=k+1call trid(a,b,c,r,u,imax)do i=1,imaxt(i)=u(i)r(i)=u(i)end dowrite(2,222) (t(i),i=1,imax)! write(*,222) (t(i),i=1,imax)if(k.lt.kmax) goto 1222 format(1x,21(2x,f18.12))endsubroutine trid(a,b,c,r,u,n)parameter (nmax=100)real gam(nmax),a(n),b(n),c(n),u(n),r(n)if (b(1)==0.) pause 'b(1)=0 in trid'bet=b(1)u(1)=r(1)/betdo j=2,ngam(j)=c(j-1)/betbet=b(j)-a(j)*gam(j)if (bet==0.) pause 'bet=0 in trid'u(j)=(r(j)-a(j)*u(j-1))/betend dodo j=n-1,1,-1u(j)=u(j)-gam(j+1)*u(j+1)end doend subroutine trid================================================== 6.3椭圆方程例题6.3.1计算程序! program overrelaxation.for! overrelaxation_iteration_methods for finite diffrence! of the 2d_Laplacian difference equationparameter(imax=25,jmax=25,pi=3.1415926)dimension T(imax,jmax),a(imax,jmax),& qx(imax,jmax),qy(imax,jmax),q(imax,jmax),an(imax,jmax) real lamda,Tu,Td,Tl,Tr,kpx,kpy,dx,dyopen(2,file='overrelaxation.dat')open(3,file='q.dat')lamda=1.5;Tu=100.;Tl=75.;Tr=50.;Td=0.;dx=1.;dy=1.;kmax=9 kpx=-0.49/2./dx;kpy=-0.49/2./dydo i=1,imaxdo j=1,jmaxT(i,j)=0.0end doend doT(1,1)=0.5*(Td+Tl) ;T(1,jmax)=0.5*(Tl+Tu)T(imax,jmax)=0.5*(Tu+Tr);T(imax,1)=0.5*(Tr+Td)do j=2,jmax-1T(1,j)=Tl; T(imax,j)=Trend dodo i=2, imax-1T(i,1)=Td; T(i,jmax)=Tuend dok=0do i=2,imax-1do j=2,jmax-1a(i,j)=T(i,j)T(i,j)=0.25*(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1))T(i,j)=lamda*T(i,j)+(1.-lamda)*a(i,j)end doend doep=abs((T(2,2)-a(2,2))/T(2,2))if(k.lt.kmax) goto 5write(*,*) 'ep=',epdo j=1,jmaxwrite(2,222) (T(i,j),i=1,imax)end do222 format(1x,25(2x,f14.10))do i=2,imax-1do j=2,jmax-1qx(i,j)=kpx*(T(i+1,j)-T(i-1,j))qy(i,j)=kpy*(T(i,j+1)-T(i,j-1))q(i,j)=sqrt(qx(i,j)*qx(i,j)+qy(i,j)*qy(i,j))if(qx(i,j).gt.0.) thenan(i,j)=atan(qy(i,j)/qx(i,j))elsean(i,j)=atan(qy(i,j)/qx(i,j))+piend ifend doend dodo j=2,jmax-1do i=2,imax-1! write(3,333) dx*(i-1),dy*(j-1),an(i,j),q(i,j)write(3,333) dx*(i-1),dy*(j-1),qx(i,j),qy(i,j)end doend do333 format(1x,4(2x,f18.12))End================================================= 例题6.3.2计算程序!!!!xadi.f90!!!!!program xadiuse AVDefuse DFLibparameter(n=19)real aj(n),bj(n),cj(n)integer statusreal T(0:n,0:n),TT(0:n,0:n)real lamda,Tu,Td,Tl,Tr,dx,dy,dt,kaopen(2,file='T3.dat')Tu=100.;Tl=75.;Tr=50.;Td=0.;dx=1.;dy=1.;kmax=10;dt=1.ka=0.835; lamda=ka*dt/(dx*dx)write(*,*) 'lamda=',lamdaT = 0.T(0,:) = Tl; T(n,:) =TrT(:,0) = Td; T(:,n) =TuT(0,0)=0.5*(Td+Tl);T(0,n)=0.5*(Tl+Tu)T(n,n)=0.5*(Tu+Tr);T(n,0)=0.5*(Tr+Td)! 系数矩阵赋值aj(:)=-lamda; ai(:)=lamdabj(:)= 2.*(1.+lamda); bi(:)=2.*(1-lamda)cj(:)=-lamda; ci(:)=lamdacall faglStartWatch(T, status)! 相当于随时间演化do k=1,kmaxdo i=1,n-1do j=1,n-1r(j)=ai(i)*T(i-1,j)+bi(i)*T(i,j)+ci(i)*T(i+1,j)end dor(1)=r(1)-aj(1)*T(i,0)r(n-1)=r(n-1)-cj(n-1)*T(i,n)call tridag(aj,bj,cj,r,u,n-1)do j=1,n-1T(i,j)=u(j)end doend do! 将矩阵T转置计算x方向call reverse(n,T,TT)do i=1,n-1do j=1,n-1r(j)=ai(i)*TT(i-1,j)+bi(i)*TT(i,j)+ci(i)*TT(i+1,j) end dor(1)=r(1)-aj(1)*TT(i,0)r(n-1)=r(n-1)-cj(n-1)*TT(i,n)call tridag(aj,bj,cj,r,u,n-1)do j=1,n-1TT(i,j)=u(j)end doend docall reverse(n,TT,T)call faglUpdate(T, status)call faglShow(T, status)end dopausecall faglClose(T, status)call faglEndWatch(T, status)do j=0,nwrite(2,222) (T(i,j),i=0,n)end do222 format(1x,20(2x,f14.10))endSUBROUTINE tridag(a,b,c,r,u,n)INTEGER n,NMAXREAL a(n),b(n),c(n),r(n),u(n)PARAMETER (NMAX=500)INTEGER jREAL bet,gam(NMAX)if(b(1).eq.0.)pause 'tridag: rewrite equations'bet=b(1)u(1)=r(1)/betdo 11 j=2,ngam(j)=c(j-1)/betbet=b(j)-a(j)*gam(j)if(bet.eq.0.)pause 'tridag failed'u(j)=(r(j)-a(j)*u(j-1))/bet11 continuedo 12 j=n-1,1,-1u(j)=u(j)-gam(j+1)*u(j+1)12 continuereturnENDSUBROUTINE REVERSE(n,T,TT) integer n,i,jreal T(0:n,0:n),TT(0:n,0:n)do i=0,ndo j=0,nTT(i,j)=T(j,i)end doend doreturnend================================================6.4 非线性偏微分方程% gasdynamincs_lw_1.mclc; clear all;xmin=-1.; xmax=1.; nx=101; h=(xmax-xmin)/(nx-1);r=0.9; cad=0.; time=0.; gamma=1.4; bcl=1.; bcr=1.;fprintf(' Space step - %f \n',h); fprintf(' \n');dd(1:nx)=zeros(1,nx); ud(1:nx)=zeros(1,nx); pd(1:nx)=zeros(1,nx);x(1:nx)=xmin+((1:nx)-1)*h;% Test 1te=' ( Test 1 )';rol=1.; ul=0.; pl=1.;ror=0.125; ur=0.; pr=0.1;tp=5.5;% Test 3%te=' ( Test 3 )';%rol=0.1; ul=0.; pl=1.;%ror=1.; ur=0.; pr=100.;%tp=0.05;% Initial conditionsdd(1:(nx-1)/2+1)=rol; ud(1:(nx-1)/2+1)=ul; pd(1:(nx-1)/2+1)=pl;dd((nx-1)/2+2:nx)=ror;ud((nx-1)/2+2:nx)=ur; pd((nx-1)/2+2:nx)=pr;[d,u,p,e,time]=lw_gasdynamics(dd,ud,pd,h,nx,r,cad,tp,gamma,bcl,bcr);% Exact solutionfor n=1:nx[density, velocity, pressure]=riemann(x(n),time,rol,ul,pl,ror,ur,pr,gamma);de(n)=density; ue(n)=velocity; pe(n)=pressure;ee(n)=pressure/((gamma-1.)*density);endfigure(1);plot(x,de,'k-',x,d,'k--');xlabel(' x ');ylabel(' density ');title([' Example 6.1 ',te]);figure(2);plot(x,ue,'k-',x,u,'k--');xlabel(' x ');ylabel(' velocity ');title([' Example 6.1 ',te]);figure(3);plot(x,pe,'k-',x,p,'k--');xlabel(' x ');ylabel(' pressure ');title([' Example 6.1 ',te]);figure(4);plot(x,ee,'k-',x,e,'k--');xlabel(' x ');ylabel(' internal energy ');title([' Example 6.1 ',te]);clear all;% Lax-Wendroff scheme (gasdynamics)function [d,u,p,e,time]=lw_gasdynamics(dd,ud,pd,h,nx,r,cad,tp,gamma,bcl,bcr)sd(1,1:nx)=dd(1:nx);sd(2,1:nx)=dd(1:nx).*ud(1:nx);sd(3,1:nx)=pd(1:nx)/(gamma-1.)+0.5*(dd(1:nx).*ud(1:nx)).*ud(1:nx);si(1:3,1:nx+1)=zeros(3,nx+1); av(1:3,1:nx)=zeros(3,nx);time=0.;while time <= tpsta=max(abs(sd(2,:)./sd(1,:))+...sqrt(gamma*(gamma-1)*(sd(3,:)./sd(1,:)-0.5*(sd(2,:).*sd(2,:))./(sd(1,:).*sd(1,:)))));tau=r*h*sqrt(1.-2.*cad)/sta;time=time+tau; %fprintf(' time - %f \n',time);% first stepif bcl == 0sl1=sd(1,2); sl2=-sd(2,2); sl3=sd(3,2);elsesl1=sd(1,2); sl2=sd(2,2); sl3=sd(3,2);endgl(1)=sl2;gl(2)=(gamma-1.)*sl3+0.5*(3.-gamma)*sl2^2/sl1;gl(3)=sl2*(gamma*sl3-0.5*(gamma-1.)*sl2^2/sl1)/sl1;for n=1:nxif (n == 1)|(n == nx)av(:,n)=[ 0.; 0.; 0.;];elseav(:,n)=cad*(sd(:,n+1)-2.*sd(:,n)+sd(:,n-1));endsr1=sd(1,n); sr2=sd(2,n); sr3=sd(3,n);gr(1)=sr2;gr(2)=(gamma-1.)*sr3+0.5*(3.-gamma)*sr2^2/sr1;gr(3)=sr2*(gamma*sr3-0.5*(gamma-1.)*sr2^2/sr1)/sr1;si(1,n)=0.5*(sl1+sr1)-(gr(1)-gl(1))*tau*0.5/h;si(2,n)=0.5*(sl2+sr2)-(gr(2)-gl(2))*tau*0.5/h;si(3,n)=0.5*(sl3+sr3)-(gr(3)-gl(3))*tau*0.5/h;gl(:)=gr(:);sl1=sr1; sl2=sr2; sl3=sr3;endif bcr == 0sr1=sd(1,nx-1); sr2=-sd(2,nx-1); sr3=sd(3,nx-1);elsesr1=sd(1,nx-1); sr2=sd(2,nx-1); sr3=sd(3,nx-1);endgr(1)=sr2;gr(2)=(gamma-1.)*sr3+0.5*(3.-gamma)*sr2^2/sr1;gr(3)=sr2*(gamma*sr3-0.5*(gamma-1.)*sr2^2/sr1)/sr1;si(1,nx+1)=0.5*(sl1+sr1)-(gr(1)-gl(1))*tau*0.5/h;si(2,nx+1)=0.5*(sl2+sr2)-(gr(2)-gl(2))*tau*0.5/h;si(3,nx+1)=0.5*(sl3+sr3)-(gr(3)-gl(3))*tau*0.5/h;% second stepfl(1)=si(2,1);fl(2)=(gamma-1.)*si(3,1)+0.5*(3.-gamma)*si(2,1)^2/si(1,1);fl(3)=si(2,1)*(gamma*si(3,1)-0.5*(gamma-1.)*si(2,1)^2/si(1,1))/si(1,1);for n=2:nx+1fr(1)=si(2,n);fr(2)=(gamma-1.)*si(3,n)+0.5*(3.-gamma)*si(2,n)^2/si(1,n);fr(3)=si(2,n)*(gamma*si(3,n)-0.5*(gamma-1.)*si(2,n)^2/si(1,n))/si(1,n);su(:,n-1)=sd(:,n-1)-(fr(:)-fl(:))*tau/h+av(:,n-1);fl(:)=fr(:);endsd(:,1:nx)=su(:,1:nx);endd(1:nx)=su(1,1:nx); u(1:nx)=su(2,1:nx)./su(1,1:nx);p(1:nx)=(gamma-1.)*(su(3,1:nx)-0.5*(su(2,1:nx).*su(2,1:nx))./su(1,1:nx));e(1:nx)=(su(3,1:nx)-0.5*(su(2,1:nx).*su(2,1:nx))./su(1,1:nx))./su(1,1:nx);return;% Solves Riemann problem for ideal gas% r1, u1, p1 initial parameters for x<0% r2, u2, p2 initial parameters for x>0% rf, uf, pf solution at the point (x,t)% c1, c2 sound speeds for x<0 and for x>0% s1, sr velocties of the fronts of S(hock)W(ave)s% shl, shr velocties of the fronts of the heads of R(arefaction)Ws% stl, str velocties of the tails of RWsfunction [rf,uf,pf]=riemann(x,t,r1,u1,p1,r2,u2,p2,gamma)g(1)=0.5*(gamma-1.)/gamma; g(2)=0.5*(gamma+1.)/gamma; g(3)=2.*gamma/(gamma-1.); g(4)=2./(gamma-1.); g(5)=2./(gamma+1.); g(6)=(gamma-1.)/(gamma+1.);g(7)=0.5*(gamma-1.); g(8)=1./gamma; g(9)=gamma-1.;s=x/t;c1 = sqrt(p1/(r1*g(8)));c2 = sqrt(p2/(r2*g(8)));du=u2-u1;% Vacuum solutionuvac=g(4)*(c1+c2)-du;if uvac <= 0.disp(' Vacuum generated by given data');return;end% Apply Newton method% Initial guessespv=0.5*(p1+p2)-0.125*du*(r1+r2)*(c1+c2);pmin=min(p1,p2); pmax=max(p1,p2); qrat=pmax/pmin;if ( qrat <= 2. ) & ( pmin <= pv & pv <= pmax )p=max(1.0e-6,pv);elseif pv < pminpnu=c1+c2-g(7)*du;pde=c1/(p1^g(1))+c2/(p2^g(1));p=(pnu/pde)^g(3);elsegel=sqrt((g(5)/r1)/(g(6)*p1+max(1.0e-6,pv)));ger=sqrt((g(5)/r2)/(g(6)*p2+max(1.0e-6,pv)));p=(gel*p1+ger*p2-du)/(gel+ger); p=max(1.0e-6,p);endendp0=p; err=1.;while err > 1.0e-5if p <= p1prat=p/p1; fl=g(4)*c1*(prat^g(1)-1.); ffl=1./(r1*c1*prat^g(2));elsea=g(5)/r1; b=g(6)*p1; qrt=sqrt(a/(b+p));fl=(p-p1)*qrt; ffl=(1.-0.5*(p-p1)/(b+p))*qrt;endif p <= p2prat=p/p2; fr=g(4)*c2*(prat^g(1)-1.); ffr=1./(r2*c2*prat^g(2));elsea=g(5)/r2; b=g(6)*p2; qrt=sqrt(a/(b+p));fr=(p-p2)*qrt; ffr=(1.-0.5*(p-p2)/(b+p))*qrt;endp=p-(fl+fr+du)/(ffl+ffr);err=2.*abs(p-p0)/(p+p0); p0=p;endif p0 < 0.p0=1.0e-6;end% Velocity of CDu0=0.5*(u1+u2+fr-fl);if s <= u0if p0 <= p1shl=u1-c1;if s <= shlrf=r1; uf=u1; pf=p1; return;elsecml=c1*(p0/p1)^g(1); stl=u0-cml;if s > stlrf=r1*(p0/p1)^g(8); uf=u0; pf=p0; return;elseuf=g(5)*(c1+g(7)*u1+s); c=g(5)*(c1+g(7)*(u1-s));rf=r1*(c/c1)^g(4); pf=p1*(c/c1)^g(3);return;endendelsepml=p0/p1; sl=u1-c1*sqrt(g(2)*pml+g(1));if s <= slrf=r1; uf=u1; pf=p1; return;elserf=r1*(pml+g(6))/(pml*g(6)+1.); uf=u0; pf=p0; return;endendelseif p0 > p2pmr=p0/p2; sr=u2+c2*sqrt(g(2)*pmr+g(1));if s >= srrf=r2; uf=u2; pf=p2; return;elserf=r2*(pmr+g(6))/(pmr*g(6)+1.); uf=u0; pf=p0; return;endelseshr=u2+c2;if s >= shrrf=r2; uf=u2; pf=p2; return;elsecmr=c2*(p0/p2)^g(1); str=u0+cmr;if s <= strrf=r2*(p0/p2)^g(8); uf=u0; pf=p0; return;elseuf=g(5)*(-c2+g(7)*u2+s); c=g(5)*(c2-g(7)*(u2-s));rf=r2*(c/c2)^g(4); pf=p2*(c/c2)^g(3);return;endendendendreturn;%%%%%%%%% gasdynamics_Godunov.m%%%%%%%%%%%%%clc; clear all;xmin=-1.; xmax=1.; nx=101; h=(xmax-xmin)/(nx-1);r=0.9; gamma=1.4; bcl=1.; bcr=1.;fprintf('gasdynamics_godunov \n'); fprintf(' Space step - %f \n',h); fprintf(' \n'); dd(1:nx-1)=zeros(1,nx-1); ud(1:nx-1)=zeros(1,nx-1); pd(1:nx-1)=zeros(1,nx-1);x(1:nx)=xmin+((1:nx)-1)*h; xa=xmin+((1:nx-1)-0.5)*h;% Test 1te=' ( Test 1 )';rol=1.; ul=0.; pl=1.;ror=0.125; ur=0.; pr=0.1;tp=15.5;% Test 2%te=' ( Test 2 )';%rol=1.; ul=0.; pl=1000.;%ror=1.; ur=0.; pr=1.;%tp=0.024;% Test 3%te=' ( Test 3 )';%rol=0.1; ul=0.; pl=1.;%ror=1.; ur=0.; pr=100.;%tp=0.05;% Test 4%te=' ( Test 4 )';%rol=5.99924; ul=19.5975; pl=460.894;%ror=5.99242; ur=-6.19633; pr=46.0950;%tp=0.05;% Test 5%te=' ( Test 5 )';%rol=3.86; ul=-0.81; pl=10.33;%ror=1.; ur=-3.44; pr=1.;%tp=0.5;% Initial conditionsdd(1:(nx-1)/2)=rol; ud(1:(nx-1)/2)=ul; pd(1:(nx-1)/2)=pl;dd((nx-1)/2+1:nx-1)=ror;ud((nx-1)/2+1:nx-1)=ur; pd((nx-1)/2+1:nx-1)=pr;[d,u,p,e,time]=godunov_gasdynamics(dd,ud,pd,h,nx,r,tp,gamma,bcl,bcr);% Exact solutionfor n=1:nx[density, velocity, pressure]=riemann(x(n),time,rol,ul,pl,ror,ur,pr,gamma);de(n)=density; ue(n)=velocity; pe(n)=pressure;ee(n)=pressure/((gamma-1.)*density);endfigure(1);plot(x,de,'k-',xa,d,'ko');xlabel(' x ');ylabel(' density ');title([' Example 6.3 ',te]);figure(2);plot(x,ue,'k-',xa,u,'ko');xlabel(' x ');ylabel(' velocity ');figure(3);plot(x,pe,'k-',xa,p,'ko');xlabel(' x ');ylabel(' pressure ');figure(4);plot(x,ee,'k-',xa,e,'ko');xlabel(' x ');ylabel(' internal energy ');clear all;% Method of S. K. Godunovfunction [d,u,p,e,time]=godunov_gasdynamics(dd,ud,pd,h,nx,r,tp,gamma,bcl,bcr)g(1)=0.5*(gamma-1.)/gamma; g(2)=0.5*(gamma+1.)/gamma; g(3)=2.*gamma/(gamma-1.); g(4)=2./(gamma-1.); g(5)=2./(gamma+1.); g(6)=(gamma-1.)/(gamma+1.);g(7)=0.5*(gamma-1.); g(8)=1./gamma; g(9)=gamma-1.;sd(1,1:nx-1)=dd(1:nx-1);sd(2,1:nx-1)=dd(1:nx-1).*ud(1:nx-1);sd(3,1:nx-1)=pd(1:nx-1)/(gamma-1.)+0.5*(dd(1:nx-1).*ud(1:nx-1)).*ud(1:nx-1);time=0.;while time <= tpsta=max(abs(sd(2,:)./sd(1,:))+...sqrt(gamma*(gamma-1)*(sd(3,:)./sd(1,:)-0.5*(sd(2,:).*sd(2,:))./(sd(1,:).*sd(1,:)))));tau=r*h/sta;time=time+tau; %fprintf(' %f \n',time);if bcl == 0[f1,f2,f3]=flux_godunov(sd(1,1),-sd(2,1),sd(3,1),sd(1,1),sd(2,1),sd(3,1),g(:));else[f1,f2,f3]=flux_godunov(sd(1,1),sd(2,1),sd(3,1),sd(1,1),sd(2,1),sd(3,1),g(:));endfleft(1)=f1; fleft(2)=f2; fleft(3)=f3;for n=1:nx-2[f1,f2,f3]=flux_godunov(sd(1,n),sd(2,n),sd(3,n),sd(1,n+1),sd(2,n+1),sd(3,n+1),g(:));fright(1)=f1; fright(2)=f2; fright(3)=f3;su(:,n)=sd(:,n)-(fright(:)-fleft(:))*tau/h;fleft(:)=fright(:);endif bcr == 0[f1,f2,f3]=flux_godunov(sd(1,nx-1),sd(2,nx-1),sd(3,nx-1),sd(1,nx-1),-sd(2,nx-1),...sd(3,nx-1),g(:));else[f1,f2,f3]=flux_godunov(sd(1,nx-1),sd(2,nx-1),sd(3,nx-1),sd(1,nx-1),sd(2,nx-1),...sd(3,nx-1),g(:));endfright(1)=f1; fright(2)=f2; fright(3)=f3;su(:,nx-1)=sd(:,nx-1)-(fright(:)-fleft(:))*tau/h;sd(:,1:nx-1)=su(:,1:nx-1);endd(1:nx-1)=su(1,1:nx-1); u(1:nx-1)=su(2,1:nx-1)./su(1,1:nx-1);p(1:nx-1)=(gamma-1.)*(su(3,1:nx-1)-0.5*(su(2,1:nx-1).*su(2,1:nx-1))./su(1,1:nx-1));e(1:nx-1)=(su(3,1:nx-1)-0.5*(su(2,1:nx-1).*su(2,1:nx-1))./su(1,1:nx-1))./su(1,1:nx-1); return;% Solves Riemann problem for ideal gas% r1, u1, p1 initial parameters for x<0% r2, u2, p2 initial parameters for x>0% rf, uf, pf solution at the point (x,t)% c1, c2 sound speeds for x<0 and for x>0% s1, sr velocties of the fronts of S(hock)W(ave)s% shl, shr velocties of the fronts of the heads of R(arefaction)Ws% stl, str velocties of the tails of RWsfunction [rf,uf,pf]=riemann(x,t,r1,u1,p1,r2,u2,p2,gamma)g(1)=0.5*(gamma-1.)/gamma; g(2)=0.5*(gamma+1.)/gamma; g(3)=2.*gamma/(gamma-1.);g(7)=0.5*(gamma-1.); g(8)=1./gamma; g(9)=gamma-1.;s=x/t;c1 = sqrt(p1/(r1*g(8)));c2 = sqrt(p2/(r2*g(8)));du=u2-u1;% Vacuum solutionuvac=g(4)*(c1+c2)-du;if uvac <= 0.disp(' Vacuum generated by given data');return;end% Apply Newton method% Initial guessespv=0.5*(p1+p2)-0.125*du*(r1+r2)*(c1+c2);pmin=min(p1,p2); pmax=max(p1,p2); qrat=pmax/pmin;if ( qrat <= 2. ) & ( pmin <= pv & pv <= pmax )p=max(1.0e-6,pv);elseif pv < pminpnu=c1+c2-g(7)*du;pde=c1/(p1^g(1))+c2/(p2^g(1));p=(pnu/pde)^g(3);elsegel=sqrt((g(5)/r1)/(g(6)*p1+max(1.0e-6,pv)));ger=sqrt((g(5)/r2)/(g(6)*p2+max(1.0e-6,pv)));p=(gel*p1+ger*p2-du)/(gel+ger); p=max(1.0e-6,p);endendp0=p; err=1.;while err > 1.0e-5if p <= p1prat=p/p1; fl=g(4)*c1*(prat^g(1)-1.); ffl=1./(r1*c1*prat^g(2));elsea=g(5)/r1; b=g(6)*p1; qrt=sqrt(a/(b+p));fl=(p-p1)*qrt; ffl=(1.-0.5*(p-p1)/(b+p))*qrt;endif p <= p2prat=p/p2; fr=g(4)*c2*(prat^g(1)-1.); ffr=1./(r2*c2*prat^g(2));a=g(5)/r2; b=g(6)*p2; qrt=sqrt(a/(b+p));fr=(p-p2)*qrt; ffr=(1.-0.5*(p-p2)/(b+p))*qrt;endp=p-(fl+fr+du)/(ffl+ffr);err=2.*abs(p-p0)/(p+p0); p0=p;endif p0 < 0.p0=1.0e-6;end% Velocity of CDu0=0.5*(u1+u2+fr-fl);if s <= u0if p0 <= p1shl=u1-c1;if s <= shlrf=r1; uf=u1; pf=p1; return;elsecml=c1*(p0/p1)^g(1); stl=u0-cml;if s > stlrf=r1*(p0/p1)^g(8); uf=u0; pf=p0; return;elseuf=g(5)*(c1+g(7)*u1+s); c=g(5)*(c1+g(7)*(u1-s));rf=r1*(c/c1)^g(4); pf=p1*(c/c1)^g(3);return;endendelsepml=p0/p1; sl=u1-c1*sqrt(g(2)*pml+g(1));if s <= slrf=r1; uf=u1; pf=p1; return;elserf=r1*(pml+g(6))/(pml*g(6)+1.); uf=u0; pf=p0; return;endendelseif p0 > p2pmr=p0/p2; sr=u2+c2*sqrt(g(2)*pmr+g(1));if s >= srrf=r2; uf=u2; pf=p2; return;rf=r2*(pmr+g(6))/(pmr*g(6)+1.); uf=u0; pf=p0; return;endelseshr=u2+c2;if s >= shrrf=r2; uf=u2; pf=p2; return;elsecmr=c2*(p0/p2)^g(1); str=u0+cmr;if s <= strrf=r2*(p0/p2)^g(8); uf=u0; pf=p0; return;elseuf=g(5)*(-c2+g(7)*u2+s); c=g(5)*(c2-g(7)*(u2-s));rf=r2*(c/c2)^g(4); pf=p2*(c/c2)^g(3);return;endendendendreturn;% Solves Riemann problem and computes fluxes for method of S. K. Godunov% r1, u1, p1 parameters in left cell% r2, u2, p2 parameters in right cell% rf, uf, pf parameters on the interface% c1, c2 sound speeds in left and right cells% s1, sr velocties of the fronts of S(hock)W(ave)s% shl, shr velocties of the fronts of the heads of R(arefaction)Ws% stl, str velocties of the tails of RWsfunction [f1,f2,f3]=flux_godunov(sl1,sl2,sl3,sr1,sr2,sr3,g)r1=sl1; u1=sl2/r1; p1=(1./g(8)-1.)*(sl3-0.5*sl2*u1);r2=sr1; u2=sr2/r2; p2=(1./g(8)-1.)*(sr3-0.5*sr2*u2);% sound speeds in left and right cellsc1 = sqrt(p1/(r1*g(8)));c2 = sqrt(p2/(r2*g(8)));du=u2-u1;% Vacuum solution。
偏微分方程数值解法试题与答案
一．填空（1553=⨯分）1．若步长趋于零时，差分方程的截断误差0→lmR ，则差分方程的解lm U 趋近于微分方程的解lm u . 此结论_______（错或对）； 2．一阶Sobolev 空间{})(,,),()(21Ω∈''=ΩL f f f y x f H y x关于内积=1),(g f _____________________是Hilbert 空间；3．对非线性（变系数）差分格式，常用 _______系数法讨论差分格式的_______稳定性； 4．写出3x y =在区间]2,1[上的两个一阶广义导数：_________________________________， ________________________________________；5．隐式差分格式关于初值是无条件稳定的. 此结论_______（错或对）。
二．（13分）设有椭圆型方程边值问题用1.0=h 作正方形网格剖分 。
（1）用五点菱形差分格式将微分方程在内点离散化； （2）用截断误差为)(2h O 的差分法将第三边界条件离散化； （3）整理后的差分方程组为 三．（12）给定初值问题xut u ∂∂=∂∂ ， ()10,+=x x u 取时间步长1.0=τ，空间步长2.0=h 。
试合理选用一阶偏心差分格式（最简显格式）, 并以此格式求出解函数),(t x u 在2.0,2.0=-=t x 处的近似值。
1．所选用的差分格式是： 2．计算所求近似值：四．（12分）试讨论差分方程()ha h a r u u r u u k l k l k l k l ττ+-=-+=++++11,1111逼近微分方程0=∂∂+∂∂xu a t u 的截断误差阶R 。
思路一：将r 带入到原式，展开后可得格式是在点（l+1/2,k+1/2）展开的。
思路二：差分格式的用到的四个点刚好是矩形区域的四个顶点，可由此构造中心点的差分格式。
计算物理学(刘金远)课后习题答案第6章：偏微分方程数值解法
第6章：偏微分方程数值解法6.1对流方程【6.1.1】考虑边值问题, 01,0(0,)0,(1,)1(,0)t x x u au x t u t u t u x x=<<>ìï==íï=î如果取：2/7x D =，(0.5),1,2,3j x j x j =-D =，8/49t D =，k t k t=D 求出111123,,u u u 【解】采用Crank-Nicolson 方法()11111111211222k k k k k k k k j j j j j j j j u u u u u u u u t x ++++-+-+éù-=-++-+ëûD D 11111113k k k k k kj j j j j j u u u u u u +++-+-+-+-=-+由边界条件：(0,)0x u t =，取100k ku u x-=D ，10,0,1,k ku u k ==L (1,)1u t =，41ku =-1 1 0 0 - (1+2s) -s 0 0 -s (1+2s) -s 0 -s (1+2s) -s 0 s L L L L 101210 0 0 0 (1-2s) s 0 0 s (1-2s) s 0 s ( 1 k n n u u s u u u +-éùéùêúêúêúêúêúêú=êúêúêúêúêúêúêúêúêúëûëûL L L L L 01211-2s) s 0 1 1kn u u u u -éùéùêúêúêúêúêúêúêúêúêúêúêúêúêúêúêúëûëûL 由初始条件：021(72j j u x j ==-，1,2,3j =，212()t s x D ==D -1 1 0 0 0-1 3 -1 0 0 0 -1 3 -1 0 -1 3 -1 0 1012340 0 0 0 01 -1 1 0 00 1 -1 1 0 1 -1 1 1 u u u u u éùéùêúêúêúêúêúêú=êúêúêúêúêúêúëûëû00123 0 1 1u u u u éùéùêúêúêúêúêúêúêúêúêúêúêúêúëûëû000117u u ==，0237u =，0357u =1112327u u -=，111000123123337u u u u u u -+-=-+=，11100234235317u u u u u -+-=-+=114591u =125191u =，136991u =6.2抛物形方程【6.2.1】分别用下面方法求定解问题22(,0)4(1)(0,)(1,)0u u t x u x x x u t u t 춶=ﶶïï=-íï==ïïî01,0x t <<>（1）取0.2x D =，1/6l =用显式格式计算1i u ；(2)取0.2,0.01x t D =D =用隐式格式计算两个时间步。
偏微分方程教程答案
偏微分方程教程答案【篇一：偏微分方程数值解习题解答案】class=txt>3页13页【篇二：3.1 常微分方程课后答案】方程dy=x+y2通过点(0,0)的第三次近似解； dx解：取?0(x)?0 12x 002xx11152x ?2(x)?y0??[x??1(x)]dx??[x?(x2)2]dx?x2?002220x1152x)]dx ?3(x)?y0??[x?(x2?0220115181x?x?x11= x2?2201604400dy 2 求方程=x-y2通过点(1,0)的第三次近似解；dx ?1(x)?y0??(x?y0)dx??xdx?x2x解：令?0(x)?012x 002xx12212152?(x)?y?[x??(x)]dx?[x?(x)]dx?x?x201?0?02220x1152x)]dx?3(x)?y0??[x?(x2?0220115181x?x?x11 =x2?2201604400则?1(x)?y0??(x?y0)dx??xdx?x2x3 题求初值问题：?dy??x2r：x??1，y?1 ?dx??y(?1)?0的解的存在区间，并求解第二次近似解，给出在解的存在空间的误差估计；b1解：因为 m=max{x2?y2}=4 则h=min(a,)= m41则解的存在区间为x?x0=x?(?1)=x?1? 4令 ?0(x)=0 ;11?1(x)=y0+?(x2?0)dx=x3+; 33x0x?2(x) 13xx4x7111312=y0+?[x?(x?)]dx=x---+ 3942186333?12x又 ?f(x,y)?2=l ?ym*l2311则：误差估计为：?2(x)??(x)?h= 24(2?1)2dy334 题讨论方程：?y在怎样的区域中满足解的存在唯一性定理的条件， dx21并求通过点（0，0）的一切解；?f(x,y)13解：因为=y在y?0上存在且连续； ?y2?23 而y3在y???0上连续 2dy33由 ?y有：y=（x+c）2 dx2131又因为y(0)=0 所以：y=x另外 y=0也是方程的解；?3?2故方程的解为：y=?x??0x?0 x?032或 y=0;6题证明格朗瓦耳不等式：设k为非负整数，f(t)和g(t)为区间??t??上的连续非负函数，且满足不等式：tf(t)?k+?f(s)g(s)ds,??t???t则有：f(t)?kexp(?g(s)ds),??t???t证明：令r（t）=?f(s)g(s)ds,则r(t)=f(?r(t)-r(t)g(t)= f(t)g(t)- r(t)g(t)?kg(t)r(t)- r(t)g(t)?kg(t);t两边同乘以exp(-?g(s)ds)则有：?tt r(t) exp(-?g(s)ds)-r(t)g(t) exp(-?g(s)ds)??t? kg(t) exp(-?g(s)ds)?两边从?到t积分：tttr(t) exp(-?g(s)ds)?-?kg(s)dsexp(-?g(r)dr)ds???tt即 r(t) ??kg(s)ds exp(-?g(r)dr)ds?stt又 f(t) ?1?k+r(t) ?k+k?g(s)exp(-?g(r)dr)ds?sts?k(1-1+ exp(-?g(r)dr)=k exp(?g(r)dr)stt即 f(t) ?k?g(r)dr;?7题假设函数f(x,y)于（x0,y0）的领域内是y的不增函数，试证方程dy= f(x,y)满足条件y(x0)= y0的解于x? x0一侧最多只有一个解；dx证明：假设满足条件y(x0)= y0的解于x? x0一侧有两个?(x),?(x)则满足：x?(x)= y0+?f(x,?(x))dxx0x?(x)= y0+?f(x,?(x))dxx0不妨假设?(x)??(x)，则?(x)- ?(x)?0xx而?(x)- ?(x)= ?f(x,?(x))dx-?f(x,?(x))dxx0x0x=?[f(x,?(x))?f(x,?(x))dxx0又因为 f(x,y)在（x0,y0）的领域内是y的增函数，则：f(x, ?(x))-f(x, ?(x))?0x则?(x)- ?(x)= ?[f(x,?(x))?f(x,?(x))dx?0x0则?(x)- ?(x)?0所以 ?(x)- ?(x)=0，即 ?(x)= ?(x) 则原命题方程满足条件y(x0)= y0的解于x? x0一侧最多只有一个解；【篇三：同济第五版高数习题答案】试说出下列各微分方程的阶数:(1)x(y′)?2yy′+x=0;解一阶.(2)xy′?xy′+y=0;解一阶.(3)xy′′′+2y′+xy=0;解三阶.(4)(7x?6y)dx+(x+y)dy=0;解一阶.(5)解二阶.;222(6) .解一阶.2. 指出下列各题中的函数是否为所给微分方程的解:(1)xy′=2y, y=5x;解y′=10x.22 因为xy′=10x=2(5x)=2y, 所以y=5x是所给微分方程的解.(2)y′+y=0, y=3sin x?4cos x;解y′=3cos x+4sin x.因为y′+y=3cos x+4sin x+3sin x?4cos x=7sin x?cos x≠0,所以y=3sin x?4cos x不是所给微分方程的解.(3)y′′?2y′+y=0, y=xe;x2xxxx2x22 解y′=2xe+xe, y′′=2e+2xe+2xe+xe=2e+4xe+xe . 因为y′′?2y′+y=2e+4xe+xe?2(2xe+xe)+xe=2e≠0,所以y=xe不是所给微分方程的解.12122x2xx2x2xxxx2xxx2x解 , .因为=0,所以是所给微分方程的解.3. 在下列各题中, 验证所给二元方程所确定的函数为所给微分方程的解:(1)(x?2y)y′=2x?y, x?xy+y=c;解将x?xy+y=c的两边对x求导得2x?y?xy′+2y y′=0,即(x?2y)y′=2x?y,所以由x?xy+y=c所确定的函数是所给微分方程的解.(2)(xy?x)y′′+xy′+yy′?2y′=0, y=ln(xy).解将y=ln(xy)的两边对x求导得再次求导得, 即 .2222222.注意到由可得 , 所以从而(xy?x)y′′+xy′+yy′?2y′=0,即由y=ln(xy)所确定的函数是所给微分方程的解.4. 在下列各题中, 确定函数关系式中所含的参数, 使函数满足所给的初始条件:(1)x?y=c, y|=5;x=0222 ,解由y|=0得0?5=c, c=?25, 故x?y=?25.x=012222 (2)y=(c+cx)e, y|=0, y′|=1;解y′=ce+2(c+cx)e. 21222xx=02x x=02x由y|=0, y′|=1得 x=0x=0,解y′=ccos(x?c).12, 即,解之得c=1, 1, 故 , 即y=?cos x .5. 写出由下列条件确定的曲线所满足的微分方程:(1)曲线在点(x, y)处的切线的斜率等于该点横坐标的平方; 解设曲线为y=y(x), 则曲线上点(x, y)处的切线斜率为y′, 由条件y′=x, 这便是所求微分方程.(2)曲线上点p(x, y)处的法线与x轴的交点为q, 且线段pq被y轴平分.解设曲线为y=y(x), 则曲线上点p(x, y)处的法线斜率为 , 由条件第pq中点的横坐标为0, 所以q点的坐标为(?x, 0), 从而有, 即yy′+2x=0.6. 用微分方程表示一物理命题: 某种气体的气压p对于温度t的变化率与气压成正比, 所温度的平方成反比.解, 其中k为比例系数.习题12?111. 试用幂级数求下列各微分方程的解:(1)y′?xy?x=1;解设方程的解为, 代入方程得,即 .可见a?1=0, 2a?a?1=0, (n+2)a1202n+2?a=0(n=1, 2, ? ? ?),n 于是,,, , , ? ? ?. , ? ? ? ,所以,即原方程的通解为 .(2)y′′+xy′+y=0;解设方程的解为, 代入方程得,即 ,于是, , ??, , , ? ? ?.所以,即原方程的通解为.(3)xy′′?(x+m)y′+my=0(m为自然数);解设方程的解为, 代入方程得,?即 .可见(a?a)m=0, (n?m)[(n+1)a01 n+1?a]=0 (n≠m),n于是a=a, , .01所以,即原方程的通解为(其中c1, c2为任意常数).(4)(1?x)y′=x2?y;解设方程的解为, 代入方程得,即可见a1+a0=0, 2a2=0, 3a3?a2?1=0, (n+1)an+1?(n?1)a n=0(n≥3),于是a1=?a0, a2=0, , (n≥4).因此原方程的通解为(c=a(5)(x+1)y′=x20为任意常数). . ?2x+y..。
(完整word版)偏微分方程数值解习题解答案
L试讨论逼近对蘇程詈+若。
的差分沁1)2)q1 二：行口匚1）解：设点为（X ? ，/曲）屮则町=讥心厶）=班勺厶+J + °（工心）（Y ）+0（F ）.ot所以截断误差为：3E=丄 ------ + ---- 「 T h 啰_喟+竺护一 o （F ）T= 0（T + 力”2）解：设点为：（X y ，/林1 ） 3则町=讥勺,_）=以E ，_+1）+ （Y ） +o （巧卩 ot “；：；=班心+1 厶+i ）=叽厶+i ）+滋（ h ）+ * 臥工心）（为 2）+o ox (X)d心;=班心亠心）=班心，/+1）+敕:;D （一力）+ 3 役;D（血 2）+0（亥2）«截断误差为：2舟A 1 ” E= ------------ + ------------ — （―+ _） T h dt dx叭：=班％厶+i ）+敗？心）（_勿+0 @2）〜dx-(史+空八dt dx 呼1_吋】+竺丛Q —O （X ）-(叱 3 +dtdx 22・试用积分插值法推导知铁。
逼近的差分裕式班勺厶叙）一班勺，乩i）+ ——-——£）dtTq2 “-” *\ | (— 4- —)dxdt = | (un t 4- un x)ds = 0* dt & \得-U] /J+U2 r+x^ A-u4 r = 0+JE (j-l? n)F (j,n)G (j^n+l)H (j-l,n+l)^% ~ 的=旳=竹“4 = W/-lMf MTh=h T-T-ll"h + LL r H + ll：4h —LL：N =Op第二章第三章第四章第五章第六章P781.如果①'（0）二0,则称工。
是』（0）的驻点（或稳定:点）-设矩阵A对称（不必正定），求证忑是』（工）的驻点.的充要条件是1心是方程加二&的解B 42・ 试用积分插值法推导知铁。
逼近的差分裕式证： 充分性：①⑻二J 缶)+ 乂(加° -b t ^+—(Ax r x)①'(Ji) = (Ax c - A, x) + A{Ax r x) aEff))S 宀沪若①0)二Q,即(山° 一氛对=0 心怎宀A X Q -h = ()目卩 Ax-b^则帀是方程Ax^b 的解卩 必要性*若心是芳程A^ = b^\解则 Ax a —h - 0 (J 4X 0 — Z?,x) = 0+^◎ (0)=(吐命-b t x) - 0+J所以町是』0)的驻点dpg%3：证明非齐次两点边值间题心現(&)二 e it (E)二 Qu与T 7面的变分间题等价：求血EH 】，认@) = G 使 J(w t ) = min J(y)其中心SiuHU (2)-d』(#) =壬仗站)-(7» —芒⑹戲(D) +而久込叭如(2.13)(提示；先把边值条件齐衩化)+d dxO 字)+梓二/ ax13页证明：令 = w(x) + v(x)其中 w(x) = Q + (x-a)0 w(a) = a yv @) = “v(a) = 0 v(^>) = 0®所以2S = 瞥+qu = j DX DX Pd r /w 血、《, 乂 、 f"丁〔P(T + ：F)]+Q(W + V )" ax dx ax* 丫 d z dv. 产 / d dw 、 豪 令 = - — O —) +(?v = /-(- —^> — +^w) = y;^ ax ax dx ax 所以(1)的等价的形式2厶” =一?0 字)= 卩ax axu(a) = a u\b) = 0a其中久=/-(-£■去字+0W )"ax ax 则由定理22知，讥是辺值间题（2）的解的充要条件是 且满定变分方程"ogf)-C/i 小 0 Vve^Pr (Zv> 一 /j )tdx + p @»： (b)f @) ① W = J(u) = J(u.+^)^— a （u^ + 兔，以.+ 无）一（/,功・ +加）［以・（E ）+加@）］ 2 □2=J(认)+ N[a@・,f)-(/,£)-+乙agd-Qfm 沁卜• Q dx dx 「（加•一/）加x +卩@加：（砂@）-卩@）戊@） Ja（3） => （4）所以可证得• 3必要性：若如 是边值间题（1）的解。
偏微分方程数值解期末试题及参考答案
《偏微分方程数值解》试卷参考答案与评分标准专业班级信息与计算科学开课系室考试日期 2006.4.14命题教师王子亭偏微分方程数值解试题(06A)参考答案与评分标准信息与计算科学专业一（10分）、设矩阵A 对称正定，定义)(),(),(21)(n R x x b x Ax x J ∈-=，证明下列两个问题等价：(1)求n R x ∈0使 )(min )(0x J x J nRx ∈=;(2)求下列方程组的解：b Ax =解: 设n R x ∈0是)(x J 的最小值点,对于任意的n R x ∈,令),(2),()()()(2000x Ax x b Ax x J x x J λλλλϕ+-+=+=, (3分)因此0=λ是)(λϕ的极小值点,0)0('=ϕ,即对于任意的n R x ∈,0),(0=-x b Ax ,特别取b Ax x -=0,则有0||||),(2000=-=--b Ax b Ax b Ax ,得到b Ax =0. (3分) 反之,若nR x ∈0满足bAx =0,则对于任意的x ,)(),(21)0()1()(00x J x Ax x x J >+==+ϕϕ,因此0x 是)(x J 的最小值点. (4分)评分标准:)(λϕ的表示式3分, 每问3分,推理逻辑性1分二（10分）、 对于两点边值问题：⎪⎩⎪⎨⎧==∈=+-=0)(,0)(),()(b u a u b a x f qu dxdu p dx d Lu 其中]),([,0]),,([,0)(min )(]),,([0min ],[1b a H f q b a C q p x p x p b a C p b a x ∈≥∈>=≥∈∈建立与上述两点边值问题等价的变分问题的两种形式：求泛函极小的Ritz 形式和Galerkin 形式的变分方程。
解: 设}0)()(),,(|{110==∈=b u a u b a H u u H 为求解函数空间,检验函数空间.取),(10b a H v ∈,乘方程两端,积分应用分部积分得到 (3分))().(),(v f fvdx dx quv dxdv dx du pv u a b a ba ==+=⎰⎰,),(10b a H v ∈∀ 即变分问题的Galerkin 形式. (3分)令⎰-+=-=b a dx fu qu dxdup u f u u a u J ])([21),(),(21)(22,则变分问题的Ritz 形式为求),(10*b a H u ∈,使)(min )(1*u J u J H u ∈= (4分)评分标准:空间描述与积分步骤3分,变分方程3分,极小函数及其变分问题4分,三（20分）、对于边值问题⎪⎩⎪⎨⎧=⨯=∈-=∂∂+∂∂∂0|)1,0()1,0(),(,12222G u G y x yux u （1）建立该边值问题的五点差分格式（五点棱形格式又称正五点格式），推导截断误差的阶。
偏微分方程数值习题解答
偏微分方程数值习题解答对于任意的)(0I C ∞∈ϕ,有⎰⎰-=bak kba dx x x f dx x x g )()()1()()()(ϕϕ则称)(x f 有k 阶广义导数,)(x g 称为)(x f 的k 阶广义导数,并记kk dxfd x g =)(注:高阶广义导数不是通过递推定义的,可能有高阶导数而没有低阶导数.2.利用)(2I L 的完全性证明))()((1I H I H m 是Hilbert 空间.证明:只证)(1I H 的完全性.设}{n f 为)(1I H 的基本列,即0||||||||||||0''01→-+-=-m n m n m n f f f f f f因此知}{},{'n n f f 都是)(2I L 中的基本列(按)(2I L 的范数).由)(2I L 的完全性,存在)(,2I L g f ∈,使 0||||,0||||0'0→-→-g f f f n n ,以下证明0||||1→-f f n (关键证明dxdfg =)由Schwarz 不等式,有00||||.|||||)())()((|ϕϕf f x x f x f n ba n -≤-⎰ 00'''|||||||||)())()((|ϕϕf f dx x x g x f nb a n -≤-⎰对于任意的)()(0I C x ∞∈ϕ,成立⎰⎰=∞→ba b a n n dx x x f dx x x f )()()()(lim ϕϕ⎰⎰=∞→ba b a nn dx x x g dx x x f )()()()(lim 'ϕϕ由⎰⎰-=ba nb a ndx x x f dx x x f )()()()(''ϕϕ取极限得到dx x x f dx x x g ba b a ⎰⎰-=)()()()('ϕϕ 即')(f x g =,即)(1I H f ∈,且0||||||||||||0''01→-+-=-f f f f f f n n n故)(1I H 中的基本列是收敛的,)(1I H 是完全的.3.证明非齐次两点边值问题证明:边界条件齐次化令)()(0a x x u -+=βα,则0u u w -=满足齐次边界条件.w 满足的方程为00Lu f Lu Lu Lw -=-=,即w 对应的边值问题为⎩⎨⎧==-=0)(,0)('b w a w Lu f Lw (P) 由定理知,问题P 与下列变分问题等价求)(min )(,**12*1w J w J H C w EHw E ∈=∈ 其中),(),(21)(0*w Lu f w w a w J --=.而Cu u a u Lu u J u u Lu f u u u u a w J +-+=-----=),(),()(~),(),(21)(000000*而200)()(),(),(C b u b p u u a u Lu +-=-β从而**)()()(~)(C b u b p u Jw J +-=β 则关于w 的变分问题P 等价于:求α=∈)(,12*a u H C u 使得)(min )()(*1u J u J a u H u α=∈=其中)()(),(),(21)(b u b p u f u u a u J β--=4就边值问题（1.2.28）建立虚功原理 解:令)(0a x u -+=βα,0u u w -=,则w 满足)(,0)('00==-=-=b w a w Lu f Lu Lu Lw等价于:1E H v ∈∀0),(),(0=--v Lu f v Lw应用分部积分,⎰⎰+-=-=-b a b a b a dx dxdvdx dw p v dx dw p vdx dx du p dx d v dx dw p dx d |)()),((还原u ,)()(),(),(),(),(),(),(),(),(000b v b p v f v u a v u a v Lu v f v u a v Lu f v w a β--=-+-=--于是,边值问题等价于:求α=∈)(,1a u H u ,使得1E H v ∈∀,成立0)()(),(),(=--b v b p v f v u a β注:形式上与用v 去乘方程两端,应用分部积分得到的相同.5试建立与边值问题等价的变分问题.解:取解函数空间为)(20I H ,对于任意)(20I H v ∈ 用v 乘方程两端,应用分部积分,得到0),(),(44=-+=-v f u dx ud v f Lu而⎰⎰-==b a b a b a dx dxdvdx u d v dx u d vdx dx u d v dx u d .|),(33334444dx dxv d dx u d dx dx vd dx u d dx dv dx u d b a b a b a ⎰⎰=+-=2222222222| 上式为),(][2222v f dx uv dxvd dx u d b a =+⎰定义dx uv dxvd dx u d v u a ba ][),(2222+=⎰,为双线性形式.变分问题为:求)(20I H u ∈,)(20I H v ∈∀),(),(v f v u a =1-41.用Galerkin Ritz -方法求边值问题⎩⎨⎧==<<=+-1)1(,0)0(102"u u x x u u 的第n 次近似)(x u n ,基函数n i x i x i ,...,2,1),sin()(==πϕ解:(1)边界条件齐次化:令x u =0,0u u w -=,则w 满足齐次边界条件,且)1(,0)0(20==-=-=w w x x Lu Lu Lw第n 次近似n w 取为∑==n i i i n c w 1ϕ,其中),...2,1(n i c i =满足的Galerkin Ritz -方程为n j x x c a j ni i j i ,...,2,1),(),(21=-=∑=ϕϕϕ 又xd jx ix ij dx x j x i dxx j x i ij dx a j i jij i ⎰⎰⎰⎰-=+=+=ππππππππϕϕϕϕϕϕ)cos()cos(2)sin()sin()cos()cos()(),(1010210''⎰-+πππjx ix sin sin 21 由三角函数的正交性,得到⎪⎩⎪⎨⎧≠=+=j i j i i a j i ,0,212),(22πϕϕ而]1)1[()(2)sin()1(),(3102--=-=-⎰jj j dx x j x x x x ππϕ 于是得到⎪⎩⎪⎨⎧+-=-=为偶数为奇数j j j j a x x c j j j j 0)1()(8),(),(2232ππϕϕϕ最后得到∑+=-+---+=]21[1233])12(1[)12(])12sin[(8)(n k n k k x k x x u ππ 2.在题1中,用0)1(=u 代替右边值条件,)(x u n 是用Galerkin Ritz -方法求解相应问题的第n 次近似,证明)(x u n 按)1,0(2L 收敛到)(x u ,并估计误差.证明:n u 对应的级数绝对收敛,由}{sin x i π的完全性知极限就是解)(x u ,其误差估计为338nR n π≤3.就边值问题(1.2.28)和基函数),...,2,1()()(n i a x x ii =-=ϕ,写出Galerkin Ritz - 方程解:边界条件齐次化,取)(0a x u -+=βα,0u u w -=, w 对应的微分方程为)(,0)('00==-=-=b w a w Lu f Lu Lu Lw对应的变分方程为 0),(),(0=--v Lu f v w a)]([)(000a x q dx dpqu dx du p dx d Lu -++-=+-=βαβ⎰⎰+-=-ba b a dx x pv b v b p v dxdp )()()(' 变分方程为dx v qu x pv b v b p v f v w a ba ⎰--+=])([)()(),(),(0'ββ取n i a x x ii,...,2,1,)()(=-=ϕ,则Galerkin -Ritz 方程为⎰⎰∑-++--+=-=ba i ba i i nj j jidxa x x q dx a x i x pb b p fc a )]()[()()()()(),(),(11βαβϕβϕϕϕ⎰+=ba j i j i j i dx q p a ][),(''ϕϕϕϕϕϕ取1,0,1===f q p ,具体计算1=n , )(1),(11a b dx a ba -==⎰ϕϕ221)(21)()()(21a b a b a b a b d -=---+-=ββ,)(211a b c -=,即解)(2101a x u u -+= 2=n :22111)()(2),(),(),(a b dx a x a a b a ba -=-=-=⎰ϕϕϕϕ3222)(34)(4),(a b dx a x a ba -=-=⎰ϕϕ3223222)(31)()()(31)(2)()(a b a b a b a b dxa x ab dx a x d ba b a -=---+-=---+-=⎰⎰ββββ 得到方程组为⎪⎪⎪⎪⎭⎫ ⎝⎛--=⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫⎝⎛----3221322)(31)(21c )(34)()(a b a b c a b a b a b a b 特别取1,0==b a ,有⎪⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛31213411121c c求解得到1,21,6131122=-=-=c c c其解为202)(21)(a x a x u u ---+=C h2 椭圆与抛物型方程有限元法§1.1 用线性元求下列边值问题的数值解:10,2sin242"<<=+-x x y y ππ0)1(,0)0('==y y此题改为4/1,0)1()0(,1"====+-h y y y y 解: 取2/1=h ,)2,1,0(==j jh x j ,21,y y 为未知数.Galerkin 形式的变分方程为),(),(v f v Lu =,其中⎰⎰+-=10210"4),(uvdx vdx u v Lu π,⎰=1)(2sin 2),(dx x xv v f π又dx v u dx v u v u vdx u ⎰⎰⎰=+-=-10''10''10'10"|因此dx uv v u v u a )4(),(12''⎰+=π在单元],[1i i i x x I -=中,应用仿射变换(局部坐标)hx x i 1--=ξ节点基函数为)3,2,1(,0,,,1)(111=⎪⎪⎪⎩⎪⎪⎪⎨⎧≤≤-=≤≤-=-=--+i other x x x h x x x x x h x x x i i i i i i i ξξξξϕ⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡-+++=++=⎰⎰⎰⎰1022210222222'111)1(41]41[]4[),(1021ξξπξξπϕπϕϕϕd h d hh dxa x x x x取2/1=h ,则计算得124),(211πϕϕ+=a122)1(41[),(210221πξξξπϕϕ+-=-+-=⎰d h h a⎰⎰-+++=10101)1)(2121(2sin )0(2sin [2),(ξξξπξξξπϕd d h h f ⎰⎰-++=1010)1(4)1(sin 2sin ξξξπξξξπd d hξξξπϕd h f ⎰+=102)2121(2sin 2),(代数方程组为⎪⎪⎭⎫ ⎝⎛=⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛),(),(),(),(),(),(212122212111ϕϕϕϕϕϕϕϕϕϕf f y y a a a a 代如求值.取4/1=h ,未知节点值为4321,,,u u u u ,方程为4,3,2,1),(),(41==∑=j f ua j i ijiϕϕϕ应用局部坐标ξ表示,⎰⎰-+++=10221022])1(41[)41(),(ξξπξξπϕϕd hh d h h a j j248]88[21022πξξπ+=+=⎰dξξξπϕϕd hh a j j ])1(41[),(1021⎰-+-=++964)1(164212πξξξπ+-=-+-=⎰d 964),(21πϕϕ+-=-j j a系数矩阵为}964,248,964{222πππ+-++-=diag A取1=f ,41)1(),(1010=-+=⎰⎰ξξξξϕd h d h f j⎰⎰-+++=+10110)1)]((2sin[2)](2sin[2),(ξξξπξξξπϕd h x h d h x h f j j j ⎰⎰-++++=1010)1)](441(2sin[21)]44(2sin[42ξξξπξξξπd j d j ⎰⎰++⨯=+++++-+=100110|)]8)1([cos(821]8)1(sin[21]8)1(sin[]8)(sin[21ξππξξπξξξπξπj d j d j j+2.就非齐次第三边值条件22'11')()(,)()(βαβα=+=+b u b u a u a u导出有限元方程.解:设方程为f qu pu Lu =+-='')( 则由),()]()[()()]()[()(),(|),)((''1122'''''v pu a u a v a p b u b v b p v pu v pu v pu b a----=-=αβαβ变分形式为:),(1b a H v ∈∀)()()()(),()()()()()()(),(),(1212''a v a p b v b p v f a v a u a p b v b u b p v qu v pu ββαα-+=-++)(),(0b u u a u u N ==记)()()()(),()()()()()()()(),(),(),(1212''a v a p b v b p v f v F a v a u a p b v b u b p v qu v pu v u A ββαα-+=-++=则上述变分形式可表示为)(),(v F v u A = 设节点基函数为),...,2,1,0)((N j x j =ϕ 则有限元方程为),...,1,0()(),(0N j F u A j Ni i j i ==∑=ϕϕϕ具体计算使用标准坐标ξ.。
偏微分方程数值习题解答
李微分方程数值解习题解答 1-1 如果0)0('=ϕ，则称0x 是)(x J 的驻点（或稳定点）.矩阵A 对称（不必正定），求证0x 是)(x J 的驻点的充要条件是：0x 是方程组 b Ax =的解证明:由)(λϕ的定义与内积的性线性性质,得),()),((21)()(0000x x b x x x x A x x J λλλλλϕ+-++=+=),(2),()(200x Ax x b Ax x J λλ+-+=),(),()(0'x Ax x b Ax λλϕ+-=必要性:由0)0('=ϕ,得,对于任何n R x ∈,有0),(0=-x b Ax ,由线性代数结论知,b Ax b Ax ==-00,0充分性: 由b Ax =0,对于任何n R x ∈,0|),(),()0(00'=+-==λλϕx Ax x b Ax即0x 是)(x J 的驻点. §1-2补充: 证明)(x f 的不同的广义导数几乎处处相等.证明:设)(2I L f ∈,)(,221I L g g ∈为)(x f 的广义导数,由广义导数的定义可知,对于任意)()(0I C x ∞∈ϕ,有⎰⎰-=ba ba dx x x f dx x x g )()()()('1ϕϕ ⎰⎰-=ba ba dx x x f dx x x g )()()()('2ϕϕ 两式相减,得到)(0)()(021I C x g g ba ∞∈∀=-⎰ϕϕ 由变分基本引理,21g g -几乎处处为零,即21,g g 几乎处处相等.补充:证明),(v u a 的连续性条件证明: 设'|)(|,|)(|M x q M x p ≤≤,由Schwarz 不等式||||.||||||||.|||||)(||),(|'''''v u M v u M dx quv v pu v u a ba +≤+=⎰11*||||.||||2v u M ≤,其中},max{'*M M M =习题:1 设)('x f 为)(x f 的一阶广义导数,试用类似的方法定义)(x f 的k 阶导数,...2,1(=k ) 解:一阶广义导数的定义,主要是从经典导数经过分部积分得到的关系式来定义,因此可得到如下定义:对于)()(2I L x f ∈,若有)()(2I L x g ∈,使得对于任意的)(0I C ∞∈ϕ,有 ⎰⎰-=bak kba dx x x f dx x x g )()()1()()()(ϕϕ则称)(x f 有k 阶广义导数,)(x g 称为)(x f 的k 阶广义导数,并记kk dxfd x g =)(注:高阶广义导数不是通过递推定义的,可能有高阶导数而没有低阶导数.2.利用)(2I L 的完全性证明))()((1I H I H m 是Hilbert 空间.证明:只证)(1I H 的完全性.设}{n f 为)(1I H 的基本列,即0||||||||||||0''01→-+-=-m n m n m n f f f f f f因此知}{},{'n n f f 都是)(2I L 中的基本列(按)(2I L 的范数).由)(2I L 的完全性,存在)(,2I L g f ∈,使0||||,0||||0'0→-→-g f f f n n ,以下证明0||||1→-f f n (关键证明dxdfg =)由Schwarz 不等式,有00||||.|||||)())()((|ϕϕf f x x f x f n ba n -≤-⎰00'''|||||||||)())()((|ϕϕf f dx x x g x f n ba n -≤-⎰对于任意的)()(0I C x ∞∈ϕ,成立⎰⎰=∞→ba ba n n dx x x f dx x x f )()()()(lim ϕϕ⎰⎰=∞→ba b a nn dx x x g dx x x f )()()()(lim 'ϕϕ由⎰⎰-=b a nba ndxxxfdxxxf)()()()(''ϕϕ取极限得到dxxxfdxxxg baba⎰⎰-=)()()()('ϕϕ即')(fxg=,即)(1IHf∈,且||||||||||||''1→-+-=-ffffffnnn故)(1IH中的基本列是收敛的,)(1IH是完全的.3.证明非齐次两点边值问题证明:边界条件齐次化令)()(axxu-+=βα,则0uuw-=满足齐次边界条件.w满足的方程为LufLuLuLw-=-=,即w对应的边值问题为⎩⎨⎧==-=0)(,0)('b w a w Lu f Lw (P) 由定理知,问题P 与下列变分问题等价求)(min )(,**12*1w J w J H C w EHw E ∈=∈ 其中),(),(21)(0*w Lu f w w a w J --=.而Cu u a u Lu u J u u Lu f u u u u a w J +-+=-----=),(),()(~),(),(21)(000000*而200)()(),(),(C b u b p u u a u Lu +-=-β从而**)()()(~)(C b u b p u Jw J +-=β 则关于w 的变分问题P 等价于:求α=∈)(,12*a u H C u使得)(min )()(*1u J u J a u H u α=∈=其中)()(),(),(21)(b u b p u f u u a u J β--=4就边值问题（）建立虚功原理 解:令)(0a x u -+=βα,0u u w -=,则w 满足)(,0)('00==-=-=b w a w Lu f Lu Lu Lw等价于:1E H v ∈∀0),(),(0=--v Lu f v Lw应用分部积分,⎰⎰+-=-=-b a b a b a dx dxdv dx dw p v dx dw p vdx dx du p dx d v dx dw p dx d |)()),(( 还原u ,)()(),(),(),(),(),(),(),(),(000b v b p v f v u a v u a v Lu v f v u a v Lu f v w a β--=-+-=--于是,边值问题等价于:求α=∈)(,1a u H u ,使得1E H v ∈∀,成立0)()(),(),(=--b v b p v f v u a β注:形式上与用v 去乘方程两端,应用分部积分得到的相同. 5试建立与边值问题等价的变分问题.解:取解函数空间为)(20I H ,对于任意)(20I H v ∈ 用v 乘方程两端,应用分部积分,得到0),(),(44=-+=-v f u dx ud v f Lu而⎰⎰-==b a b a b a dx dxdvdx u d v dx u d vdx dx u d v dx u d .|),(33334444 dx dxv d dx u d dx dx vd dx u d dx dv dx u d b a b a b a ⎰⎰=+-=2222222222| 上式为),(][2222v f dx uv dxvd dx u d b a =+⎰定义dx uv dxvd dx u d v u a ba ][),(2222+=⎰,为双线性形式.变分问题为:求)(20I H u ∈,)(20I H v ∈∀),(),(v f v u a =1-41.用Galerkin Ritz -方法求边值问题⎩⎨⎧==<<=+-1)1(,0)0(102"u u x x u u 的第n 次近似)(x u n ,基函数n i x i x i ,...,2,1),sin()(==πϕ解:(1)边界条件齐次化:令x u =0,0u u w -=,则w 满足齐次边界条件,且)1(,0)0(20==-=-=w w x x Lu Lu Lw第n 次近似n w 取为∑==n i i i n c w 1ϕ,其中),...2,1(n i c i =满足的Galerkin Ritz -方程为n j x x c a j ni i j i ,...,2,1),(),(21=-=∑=ϕϕϕ 又xd jx ix ij dx x j x i dxx j x i ij dx a j i jij i ⎰⎰⎰⎰-=+=+=ππππππππϕϕϕϕϕϕ)cos()cos(2)sin()sin()cos()cos()(),(1010210''⎰-+πππjx ix sin sin 21由三角函数的正交性,得到⎪⎩⎪⎨⎧≠=+=j i j i i a j i ,0,212),(22πϕϕ而]1)1[()(2)sin()1(),(3102--=-=-⎰jj j dx x j x x x x ππϕ 于是得到⎪⎩⎪⎨⎧+-=-=为偶数为奇数j j j j a x x c j j j j 0)1()(8),(),(2232ππϕϕϕ最后得到∑+=-+---+=]21[1233])12(1[)12(])12sin[(8)(n k n k k x k x x u ππ 2.在题1中,用0)1(=u 代替右边值条件,)(x u n 是用Galerkin Ritz -方法求解相应问题的第n 次近似,证明)(x u n 按)1,0(2L 收敛到)(x u ,并估计误差. 证明:n u 对应的级数绝对收敛,由}{sin x i π的完全性知极限就是解)(x u ,其误差估计为338nR n π≤3.就边值问题和基函数),...,2,1()()(n i a x x i i =-=ϕ,写出Galerkin Ritz -方程解:边界条件齐次化,取)(0a x u -+=βα,0u u w -=, w 对应的微分方程为)(,0)('00==-=-=b w a w Lu f Lu Lu Lw对应的变分方程为0),(),(0=--v Lu f v w a)]([)(000a x q dx dpqu dx du p dx d Lu -++-=+-=βαβ⎰⎰+-=-ba b a dx x pv b v b p v dxdp )()()(' 变分方程为dx v qu x pv b v b p v f v w a ba ⎰--+=])([)()(),(),(0'ββ取n i a x x i i ,...,2,1,)()(=-=ϕ,则Galerkin -Ritz 方程为⎰⎰∑-++--+=-=ba i ba i i nj j jidxa x x q dx a x i x pb b p fc a )]()[()()()()(),(),(11βαβϕβϕϕϕ⎰+=ba j i j i j i dx q p a ][),(''ϕϕϕϕϕϕ取1,0,1===f q p ,具体计算1=n , )(1),(11a b dx a ba -==⎰ϕϕ221)(21)()()(21a b a b a b a b d -=---+-=ββ,)(211a b c -=,即解)(2101a x u u -+= 2=n :22111)()(2),(),(),(a b dx a x a a b a ba -=-=-=⎰ϕϕϕϕ3222)(34)(4),(a b dx a x a ba -=-=⎰ϕϕ3223222)(31)()()(31)(2)()(a b a b a b a b dxa x ab dx a x d ba b a -=---+-=---+-=⎰⎰ββββ 得到方程组为⎪⎪⎪⎪⎭⎫ ⎝⎛--=⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫⎝⎛----3221322)(31)(21c )(34)()(a b a b c a b a b a b a b特别取1,0==b a ,有⎪⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛31213411121c c求解得到1,21,6131122=-=-=c c c其解为202)(21)(a x a x u u ---+=C h2 椭圆与抛物型方程有限元法§ 用线性元求下列边值问题的数值解:10,2sin242"<<=+-x x y y ππ0)1(,0)0('==y y此题改为4/1,0)1()0(,1"====+-h y y y y解: 取2/1=h ,)2,1,0(==j jh x j ,21,y y 为未知数.Galerkin 形式的变分方程为),(),(v f v Lu =,其中⎰⎰+-=10210"4),(uvdx vdx u v Lu π,⎰=1)(2sin 2),(dx x xv v f π又dx v u dx v u v u vdx u ⎰⎰⎰=+-=-10''10''10'10"|因此dx uv v u v u a )4(),(12''⎰+=π在单元],[1i i i x x I -=中,应用仿射变换(局部坐标)hx x i 1--=ξ节点基函数为)3,2,1(,0,,,1)(111=⎪⎪⎪⎩⎪⎪⎪⎨⎧≤≤-=≤≤-=-=--+i other x x x h x x x x x h x x x i i i i i i i ξξξξϕ⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡-+++=++=⎰⎰⎰⎰1022210222222'111)1(41]41[]4[),(1021ξξπξξπϕπϕϕϕd h d hh dxa x x x x取2/1=h ,则计算得124),(211πϕϕ+=a122)1(41[),(210221πξξξπϕϕ+-=-+-=⎰d h h a⎰⎰-+++=10101)1)(2121(2sin )0(2sin [2),(ξξξπξξξπϕd d h h f ⎰⎰-++=1010)1(4)1(sin 2sin ξξξπξξξπd d hξξξπϕd h f ⎰+=102)2121(2sin 2),(代数方程组为⎪⎪⎭⎫ ⎝⎛=⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛),(),(),(),(),(),(212122212111ϕϕϕϕϕϕϕϕϕϕf f y y a a a a 代如求值.取4/1=h ,未知节点值为4321,,,u u u u ,方程为4,3,2,1),(),(41==∑=j f ua j i ijiϕϕϕ应用局部坐标ξ表示,⎰⎰-+++=10221022])1(41[)41(),(ξξπξξπϕϕd hh d h h a j j248]88[21022πξξπ+=+=⎰dξξξπϕϕd hh a j j ])1(41[),(1021⎰-+-=++964)1(164212πξξξπ+-=-+-=⎰d 964),(21πϕϕ+-=-j j a系数矩阵为}964,248,964{222πππ+-++-=diag A取1=f ,41)1(),(1010=-+=⎰⎰ξξξξϕd h d h f j⎰⎰-+++=+10110)1)]((2sin[2)](2sin[2),(ξξξπξξξπϕd h x h d h x h f j j j ⎰⎰-++++=1010)1)](441(2sin[21)]44(2sin[42ξξξπξξξπd j d j⎰⎰++⨯=+++++-+=100110|)]8)1([cos(821]8)1(sin[21]8)1(sin[]8)(sin[21ξππξξπξξξπξπj d j d j j+2.就非齐次第三边值条件22'11')()(,)()(βαβα=+=+b u b u a u a u导出有限元方程.解:设方程为f qu pu Lu =+-='')( 则由),()]()[()()]()[()(),(|),)((''1122'''''v pu a u a v a p b u b v b p v pu v pu v pu b a----=-=αβαβ变分形式为:),(1b a H v ∈∀)()()()(),()()()()()()(),(),(1212''a v a p b v b p v f a v a u a p b v b u b p v qu v pu ββαα-+=-++)(),(0b u u a u u N ==记)()()()(),()()()()()()()(),(),(),(1212''a v a p b v b p v f v F a v a u a p b v b u b p v qu v pu v u A ββαα-+=-++=则上述变分形式可表示为)(),(v F v u A =设节点基函数为),...,2,1,0)((N j x j =ϕ 则有限元方程为),...,1,0()(),(0N j F u A j Ni i j i ==∑=ϕϕϕ具体计算使用标准坐标ξ.。
偏微分课后习题答案终极版
∂u t sΔxΔt ∂t
cρ
∂u ∂ 2u s Δ x Δ t = k t ∂t ∂x 2
x sΔxΔt
−
4k1 (u − u1 )sΔxΔt l
消去 sΔxΔt ，再令 Δx → 0 ， Δt → 0 得精确的关系：
cρ
∂ 2 u 4k ∂u = k 2 − 1 (u − u1 ) l ∂t ∂x
πl 2
4
为 S 。由假设，在任意时刻 t 到 t + Δt 内流入截面坐标为 x 到 x + Δx 一小段细杆的热量为
dQ1 = k
∂u ∂x
x + Δx sΔt − k
∂u ∂x
x sΔt = k
∂ 2u ∂x 2
x sΔxΔt
杆表面和周围介质发生热交换，可看作一个“被动”的热源。由假设，在时刻 t 到 t + Δt 在 截面为 x 到 x + Δx 一小段中产生的热量为
2
代入化简即得 (3) 因
x 2uηη = 0 uηη = 0 ( x ≠ 0)
u xx − 2 cos xu xy − (3 + sin 2 x)u yy − yu y = 0
Δ = cos 2 x + (3 + sin 2 x) = 4 > 0 为双曲型.特征方程为 dy dy ( ) 2 + 2 cos x − (3 + sin 2 x) = 0 dx dx
于是得运动方程
ρ ( x) s( x) ⋅ Δx ⋅ utt ( x, t ) = ESu x ( x + Δx) | x + Δx − ESu x ( x) | x ρ ( x) s( x)u tt =
(完整word版)偏微分方程数值解法答案
1. 课本2p 有证明2. 课本812,p p 有说明3. 课本1520,p p 有说明4. Rit2法，设n u 是u 的n 维子空间，12,...n ϕϕϕ是n u 的一组基底，n u 中的任一元素n u 可表为1nn i i i u c ϕ==∑，则,1111()(,)(,)(,)(,)22j nnn n n n i j i j j i j j J u a u u f u a c c c f ϕϕϕ===-=-∑∑是12,...n c c c 的二次函数，(,)(,)i j j i a a ϕϕϕϕ=，令()0n jJ u c ∂=∂，从而得到12,...n c c c 满足1(,)(,),1,2...niji j i a c f j n ϕϕϕ===∑，通过解线性方程组，求的i c ，代入1nn i i i u c ϕ==∑，从而得到近似解n u 的过程称为Rit2法简而言之，Rit2法：为得到偏微分方程的有穷维解，构造了一个近似解，1nn i ii u c ϕ==∑，利用,1111()(,)(,)(,)(,)22j nnn n n n i j i j j i j j J u a u u f u a c c c f ϕϕϕ===-=-∑∑确定i c ，求得近似解n u 的过程Galerkin 法：为求得1nn i ii u c ϕ==∑形式的近似解，在系数i c 使n u 关于n V u ∈，满足(,)(,)n a u V f V =，对任意nV u ∈或（取,1j V j nϕ=≤≤）1(,)(,),1,2...nijij i a cf j n ϕϕϕ===∑的情况下确定i c ，从而得到近似解1nn i i i u c ϕ==∑的过程称Galerkin 法为 Rit2-Galerkin 法方程：1(,)(,)nijij i a cf ϕϕϕ==∑5. 有限元法：将偏微分方程转化为变分形式，选定单元的形状，对求解域作剖分，进而构造基函数或单元形状函数，形成有限元空间，将偏微分方程转化成了有限元方程，利用有效的有限元方程的解法，给出偏微分方程近似解的过程称为有限元法。
经典偏微分方程课后习题答案
第四章 抛物型微分方程有限差分法1设已知初边值问题22, 01, 0<(,0)sin , 01(0,)(1,)0, 0 u ux t t x u x x x u t u t t T π⎧∂∂=<<⎪∂∂⎪⎪=≤≤⎨⎪==≤≤⎪⎪⎩T ≤, 试用最简显格式求上述问题的数值解。
取h=0.1,r=0.1.0 1/10 2/10 … 1 T 2τ τt解: 1.矩形网格剖分区域. 取空间步长1, 时间2510h =0.00τ=以及0.01τ=的矩形网格剖分区域, 用节点)表示坐标点(,j k (,)(,)j k x t jh k τ=, 0,1,...1/; 0,1,...,/j h k T τ==, 如图所示.显然, 我们需要求解这(1/1)(/1)h T τ+×+个点对应的函数值. 事实上由已知初边界条件蓝标附近的点可直接得到, 所以只要确定微分方程的解在其它点上的取值即可. 沿用记号[]k(,)j j k u x t =。
u 2. 建立差分格式, 对于11,...1; 0,1,...,1Tj k hτ=−=−, 用向前差商代替关于时间的一阶偏导数, 用二阶中心差商代替关于空间的二阶偏导数, 则可定义最简显格式:1122k k k k k1jj j j u u u u u h ++−+=. 变形j τ−−有:1112(12) (k k k kj j j j u ru r u ru r h τ+−+=+−+=(4.1)用向后差商代替关于时间的一阶偏导数, 用二阶中心差商代替关于空间的二阶偏导数, 则可定义最简显格式最简隐格式:111122k k k k k j jj j j u u u u u h τ++++−−+=11+−1kj +,变形有:1111(12) k k k j j j ru r u ru u ++−−−++−= (4.2)(4.1)*0.5+(4.2)*0.5得CN 格式为:111112222k k k k k k k k j jj j j j j j u u u u u u u u h τ+++−+−−++−+=111++−1kj +x x变形有:111111(22)(22) k k k k k j j j j j ru r u ru ru r u ru ++−−+−−++−=+−+ (4.3)3 初边界点差分格式处理.对于初始条件u x (,0)sin , 01=π≤≤h 离散为(4.4)0sin 0,1,...1/j u jh j π==对于边界条件离散为(0,)(1,)0, 0 u t u t t T ==≤≤00 0,1,.../k k N u u k T τ===(4.5)总结: 联立方程(4.1)(4.4)(4.5)得到已知问题的最简显格式差分方程组:11100(12)1 1,...1; 0,1,...,1sin 0,1,...1/0 0,1,.../k k k k j j j j jk k N u ru r u ru T j k h u jh j h u u k T τπτ+−+⎧=+−+⎪⎪=−=−⎪⎨⎪==⎪⎪===⎩ 联立方程(4.2)( 4.4)( 4.5)得到已知问题的最简隐格式差分方程组:1111100(12) 1 1,...1; 0,1,...,1sin 0,1,...1/0 0,1,.../k k k k j j j j jk k N ru r u ru u T j k h u jh j h u u k T τπτ++−−+⎧−++−=⎪⎪=−=−⎪⎨⎪==⎪⎪===⎩ 联立方程(4.3)( 4.4)( 4.5)得到已知问题的CN 格式差分方程组:11111100(22)(22) 1 1,...1; 0,1,...,1sin 0,1,...1/0 0,1,.../k k k k k j j j j j jk k N ru r u ru ru r u ru T j k h u jh j h u u k T τπτ++−−+−⎧−++−=+−+⎪⎪=−=−⎪⎨⎪==⎪⎪===⎩1k j + 4 求解并显示结果利用软件计算(Matlab)如上最简显格式差分方程组.h=1/10;tau=0.0025;T=0.5; r=tau/h^2;M=1/h+1;N=T/tau+1; u=zeros(M,N);for m=1:Mu(m,1)=sin((m-1)*h*pi); endu(1,1:N)=0;u(M,1:N)=0;for n=1:N-1for m=2:M-1u(m,n+1)=r*(u(m+1,n)+u(m-1,n))+(1-2*r)*u(m,n); end end u=u’ 这样我们就计算出不同时刻不同位置k t j x 对应的函数值(,)j k u x t 取tau=0.0025, 即r=0.25绘图, 取tau=0.01, r=1再绘图,如图()图4.2 习题1数值解图示(左r=0.25, 右r=1)2.试构造初边值问题 ()()()()(), 0.51, 0,,0, 0.51,0.5,0, 1,0.51,, 0u u x x x T t x x u x x x u ⎪∂u t t u t t T x ϕ⎧∂∂∂⎛⎞=<<<≤⎜⎟⎪∂∂∂⎝⎠⎪⎪=≤≤⎨⎪==−≤≤⎪∂⎩的显格式，并给出其按最大范数稳定的充分条件。
偏微分方程数值解法题解
偏微分方程数值解法（带程序）例1 求解初边值问题22,(0,1),012,(0,]2(,0)12(1),[,1)2(0,)(1,)0,0u ux t t x x x u x x x u t u t t ⎧⎪⎪⎪⎧⎨⎪⎪⎪⎨⎪⎪⎪⎪⎩⎩∂∂=∈>∂∂∈=-∈==>要求采用树脂格式 111(2)n n n n nj j j j j u u u u u λ++-=+-+，2()tx λ∆=∆，完成下列计算： （1） 取0.1,0.1,x λ∆==分别计算0.01,0.02,0.1,t =时刻的数值解。
（2） 取0.1,0.5,x λ∆==分别计算0.01,0.02,0.1,t =时刻的数值解。
（3） 取0.1, 1.0,x λ∆==分别计算0.01,0.02,0.1,t =时刻的数值解。
并与解析解22()22181(,)sin()sin()2n t n u n x t n x e n ππππ∞-==∑进行比较。
解：程序function A=zhongxinchafen(x,y,la) U=zeros(length(x),length(y)); for i=1:size(x,2)if x(i)>0&x(i)<=0.5 U(i,1)=2*x(i); elseif x(i)>0.5&x(i)<1 U(i,1)=2*(1-x(i)); end endfor j=1:length(y)-1for i=1:length(x)-2U(i+1,j+1)=U(i+1,j)+la*(U(i+2,j)-2*U(i+1,j)+U(i,j)); end endA=U(:,size(U,2))function u=jiexijie1(x,t) for i=1:size(x,2) k=3;a1=(1/(1^2)*sin(1*pi/2)*sin(1*pi*x(i))*exp(-1^2*pi^2*t));a2=a1+(1/(2^2)*sin(2*pi/2)*sin(2*pi*x(i))*exp(-2^2*pi^2*t));while abs(a2-a1)>0.00001a1=a2;a2=a1+(1/(k^2)*sin(k*pi/2)*sin(k*pi*x(i))*exp(-k^2*pi^2*t));k=k+1;endu(i)=8/(pi^2)*a2;endclc; %第1题第1问clear;t1=0.01;t2=0.02;t3=0.1;x=[0:0.1:1];y1=[0:0.001:t1];y2=[0:0.001:t2];y3=[0:0.001:t3];la=0.1;subplot(131)A1=zhongxinchafen(x,y1,la);u1=jiexijie1(x,t1)line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold online(x,u1,'color','b','linewidth',1);A2=zhongxinchafen(x,y2,la);u2=jiexijie1(x,t2)line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,u2,'color','b','linewidth',1);A3=zhongxinchafen(x,y3,la);u3=jiexijie1(x,t3)line(x,A3,'color','r','linestyle',':','linewidth',1.5);line(x,u3,'color','b','linewidth',1); title('例1(1)');subplot(132);line(x,u1,'color','b','linewidth',1);line(x,u2,'color','b','linewidth',1);line(x,u3,'color','b','linewidth',1);title('解析解');subplot(133);line(x,A1,'color','r','linestyle',':','linewidth',1.5);line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,A3,'color','r','linestyle',':','linewidth',1.5);title('数值解');clc; %第1题第2问clear;t1=0.01;t2=0.02;t3=0.1;x=[0:0.1:1];y1=[0:0.005:t1];y2=[0:0.005:t2];y3=[0:0.005:t3];la=0.5;subplot(131);A1=zhongxinchafen(x,y1,la);u1=jiexijie1(x,t1)line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold online(x,u1,'color','b','linewidth',1);A2=zhongxinchafen(x,y2,la);u2=jiexijie1(x,t2)line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,u2,'color','b','linewidth',1);A3=zhongxinchafen(x,y3,la);u3=jiexijie1(x,t3)line(x,A3,'color','r','linestyle',':','linewidth',1.5); line(x,u3,'color','b','linewidth',1);title('例1(2)'); subplot(132);line(x,u1,'color','b','linewidth',1); line(x,u2,'color','b','linewidth',1);line(x,u3,'color','b','linewidth',1);title('解析解'); subplot(133);line(x,A1,'color','r','linestyle',':','linewidth',1.5); line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,A3,'color','r','linestyle',':','linewidth',1.5);title('数值解');clc; %第1题第3问 clear;t1=0.01;t2=0.02;t3=0.1;x=[0:0.1:1];y1=[0:0.01:t1];y2=[0:0.01:t2];y3=[0:0.01:t3];la=1.0; subplot(131);A1=zhongxinchafen(x,y1,la);u1=jiexijie1(x,t1)line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1);A2=zhongxinchafen(x,y2,la);u2=jiexijie1(x,t2) line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,u2,'color','b','linewidth',1);A3=zhongxinchafen(x,y3,la);u3=jiexijie1(x,t3) line(x,A3,'color','r','linestyle',':','linewidth',1.5); line(x,u3,'color','b','linewidth',1);title('例1(3)'); subplot(132);line(x,u1,'color','b','linewidth',1); line(x,u2,'color','b','linewidth',1);line(x,u3,'color','b','linewidth',1);title('解析解'); subplot(133);line(x,A1,'color','r','linestyle',':','linewidth',1.5); line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,A3,'color','r','linestyle',':','linewidth',1.5);title('数值解'); 运行结果：表1：取0.1,0.1,x λ∆==0.01,0.02,0.1,t =时刻的解析解与数值解表2：取0.1,0.5,x λ∆==0.01,0.02,0.1,t =时刻的解析解与数值解表3：取0.1, 1.0,x λ∆==0.01,0.02,0.1,t =时刻的解析解与数值解图1：取0.1,0.1,x λ∆==0.01,0.02,0.1,t =时刻的解析解与数值解图2：取0.1,0.5,x λ∆==0.01,0.02,0.1,t =时刻的解析解与数值解图3：取0.1, 1.0,x λ∆==0.01,0.02,0.1,t =时刻的解析解与数值解例2 用Crank-Nicolson 格式完成例1的所有任务。
偏微分方程习题及答案
偏微分方程习题及答案【篇一：偏微分方程数值解法期末考试题答案】题答案及评分标准学年学期：专业：班级：课程：教学大纲：使用教材：教材作者：出版社：数学与应用数学数学偏微分方程数值解法《偏微分方程数值解法》教学大纲（自编，2006）《偏微分方程数值解法》陆金甫、关治清华大学出版社一、判断题（每小题1分，共10分）1、（o）2、（o）3、（x）4、（x）5、（o）6、（o）7、（o）8、（x）9、（x） 10、（o）二、选择题（每小题2分，共10分） 11、（d） 12、（a） 13、（c） 14、（b）15、（c）三、填空题（每小题2分，共20分）?2?216、2?2??x1?x2?2?2 17、a=[4 5 9;23 5 17;11 23 1] 18、y=exp(-t/3)*sin(3*t) ?xn19、help 20、zeros(m,n)21、inva(a)*b或者a/b22、a=sym([cos(x-y),sin(x+y);exp(x-y),(x-1)^3])?(s)?1?(s)?c[??(s)]2?023、a[?2(s)]2?2b?224????v(?)ed? 25、i?xu(xj,tn?1)?u(xj,tn)?四、计算题：（每小题12分，共36分）?u?u?0（x?r,t?0）的有限差分方程（两层显示26、写成对流方程?a?t?x格式，用第n层计算第n+1层），并把有限差分方程改写为便于计算的迭代格式???/h为网格比。
解：在点(xj,tn)处，差分方程为?1un?unjj??anunj?1?ujh?0（j?0,?1,?2,，n?0,1,2,）（8分）便于计算的形式为?1nnn???/h （4分） un?u?a?(u?ujjj?1j)，?u?2u?a2的有限差分方程（中心差分格式，用第n层27、写出扩散方程?t?x计算第n+1层），并把有限差分方程改写为便于计算的迭代格式，???/h2为网格比。
偏微分方程数值解例题答案
二、改进的Euler 方法梯形方法的迭代公式(1.10)比Euler 方法精度高,但其计算较复杂,在应用公式(1.10)进行计算时,每迭代一次,都要重新计算函数),(y x f 的值,且还要判断何时可以终止或转下一步计算.为了控制计算量和简化计算法,通常只迭代一次就转入下一步计算.具体地说,我们先用Euler 公式求得一个初步的近似值1+n y ,称之为预测值,然后用公式(1.10)作一次迭代得1+n y ,即将1+n y 校正一次.这样建立的预测－校正方法称为改进的Euler 方法:预测: ),,(1n n n n y x hf y y +=+ 校正:)].,(),([2111+++++=n n n n n n y x f y x f hy y(1.15)这个计算公式也可以表示为11(,),(,),1().2p n n nc n n p n p cy y hf x y y y hf x y y y y ++⎧=+⎪⎪=+⎪⎨⎪=+⎪⎪⎩例1 取步长0.1h =，分别用Euler 方法及改进的Euler 方法求解初值问题d (1),01,d (0) 1.yy xy x xy ⎧=-+≤≤⎪⎨⎪=⎩ 解 这个初值问题的准确解为()1(21)xy x e x =--. 根据题设知).1(),(xy y y x f +-=(1) Euler 方法的计算式为)],1([1.01n n n n n y x y y y +⨯-=+由1)0(0==y y , 得,9.0)]101(1[1.011=⨯+⨯⨯-=y,8019.0)]9.01.01(9.0[1.09.02=⨯+⨯⨯-=y这样继续计算下去，其结果列于表9.1.(2) 改进的Euler 方法的计算式为110.1[(1)],0.1[(1)],1(),2p n n n n c n p n p n p c y y y x y y y y x y y y y ++⎧=-⨯+⎪=-⨯+⎪⎪⎨⎪=+⎪⎪⎩由1)0(0==y y ，得110.1[1(101)]0.9,10.1[0.9(10.10.9)]0.9019,1(0.90.9019)0.900952p c y y y ⎧=-⨯⨯+⨯=⎪⎪=-⨯⨯+⨯=⎨⎪⎪=+=⎩ 20.900950.1[0.90095(10.10.90095)]0.80274,0.900950.1[0.80274(10.20.80274)]0.80779,1(0.802740.80779)0.805262p c y y y ⎧=-⨯⨯+⨯=⎪⎪=-⨯⨯+⨯=⎨⎪⎪=+=⎩ 这样继续计算下去，其结果列于表9.1.从表9.1可以看出,Euler 方法的计算结果只有2位有效数字,而改进的Euler 方法确有3位有效数字,这表明改进的Euler 方法的精度比Euler 方法高.例2 试用Euler 方法、改进的Euler 方法及四阶经典R-K 方法在不同步长下计算初值问题⎪⎩⎪⎨⎧=≤≤+-=1)0(,10),1(d d y x xy y xy在0.2、0.4、0.8、1.0处的近似值，并比较它们的数值结果.解 对上述三种方法,每执行一步所需计算)1(),(xy y y x f +-=的次数分别为1、2、4。
最新偏微分方程数值解试题参考答案
偏微分方程数值解一（10分）、设矩阵A 对称正定，定义)(),(),(21)(n R x x b x Ax x J ∈-=，证明下列两个问题等价：(1)求n R x ∈0使)(min )(0x J x J n Rx ∈=;(2)求下列方程组的解：b Ax = 解: 设n R x ∈0是)(x J 的最小值点,对于任意的n R x ∈,令),(2),()()()(2000x Ax x b Ax x J x x J λλλλϕ+-+=+=, (3分)因此0=λ是)(λϕ的极小值点,0)0('=ϕ,即对于任意的n R x ∈,0),(0=-x b Ax ,特别取b Ax x -=0,则有0||||),(2000=-=--b Ax b Ax b Ax ,得到b Ax =0. (3分) 反之,若nR x ∈0满足bAx =0,则对于任意的x ,)(),(21)0()1()(00x J x Ax x x J >+==+ϕϕ,因此0x 是)(x J 的最小值点. (4分)评分标准:)(λϕ的表示式3分, 每问3分,推理逻辑性1分二（10分）、对于两点边值问题：⎪⎩⎪⎨⎧==∈=+-=0)(,0)(),()(b u a u b a x f qu dxdu p dx d Lu 其中]),([,0]),,([,0)(min )(]),,([0min ],[1b a H f q b a C q p x p x p b a C p b a x ∈≥∈>=≥∈∈建立与上述两点边值问题等价的变分问题的两种形式：求泛函极小的Ritz 形式和Galerkin 形式的变分方程。
解: 设}0)()(),,(|{11==∈=b u a u b a H u u H 为求解函数空间,检验函数空间.取),(10b a H v ∈,乘方程两端,积分应用分部积分得到 (3分))().(),(v f fvdx dx quv dxdv dx du p v u a b a ba ==+=⎰⎰,),(1b a H v ∈∀ 即变分问题的Galerkin 形式. (3分)令⎰-+=-=b a dx fu qu dxdup u f u u a u J ])([21),(),(21)(22,则变分问题的Ritz 形式为求),(1*b a H u ∈,使)(m in )(10*u J u J H u ∈= (4分) 评分标准:空间描述与积分步骤3分,变分方程3分,极小函数及其变分问题4分,三（20分）、对于边值问题⎪⎩⎪⎨⎧=⨯=∈-=∂∂+∂∂∂0|)1,0()1,0(),(,12222G u G y x yux u （1）建立该边值问题的五点差分格式（五点棱形格式又称正五点格式），推导截断误差的阶。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
X
(1
j
分基函数的取法
其中Rij是以(Xj,yj)为顶点的矩形单元x, y0为&的底和商的长度。
8.何为三角剖分，基函数怎样取？
三角剖分：设G是多边形域(否则可用多边形域逼近它)，将G分割成有限个三角形之
和，使不同三角形无重叠的部，且任一三角形的顶点不属于其他三角形的部，这样就把G
10，11题不会。 在此将14题推导过程介绍如下：
12.对Possion方程f (x, y)，建立五点差分格式，并估计截断误差。
取定沿x轴和y轴方向的步长h1和h2,沿x,y方向分别用二阶中心差商代替，则
2 2
i1,ji,j i 1,j i,j 1 i, j i,j 1
h ij[22]fij(五点差分格式)
在Richardson格式(4.1.10)中以 :—(n1
j
）代入，便得Du Fort Fran：el格式:
；1(n1
h
n
1n
12
2
1
j
j
t
2
t2
6
n
1n
12
2
1
j
j
t
2
t2
6
n 1
n 1
①-
②得-
j
j
1
2
t
6
从而得到了微分方程左边的误差
同理可得微分方程右边的误差:
2.44
旦洁—h
h
1
3
3
t3
144
vk+1exp(ia jh) = r exp(ia( j +1)h)+(1-2r)exp( ia jh) + r exp(ia( j-1)h)vk消 去
exp(iajh)则知增长因子
2ah C1=(xp,t)=(1-2r) + r(exp(iah)+exp(-iah))=1-2r(1-cosah)=1-4rsin乙由于 霁甲)在［0，n］中分布稠密，(随h?0)为使Cgt)满足vonNeu-Mann条件，
2(x, yj)
h；
4(Xi,yj)
h：
6(x, yj)
(h6)
2
y
12
4
y
360
6
y
("b)
截断误差为
4
4
Rj()
(Xi, yj)
n(Xi, yj)
12
(h14
h；
4)ij
(h4)
(h2)
12x
y
13.对possion方程建立，极坐标形式的差分格式
poission方程的极坐标形式为
2
—]
f(
其中
tan
nn
a(i,j)Ci(f,j), j1,2...n，通过解线性方程组，求的Ci，代入unCi i，
i 1i 1
从而得到近似解Un的过程称为Rit2法
n
简而言之，Rit2法：为得到偏微分方程的有穷维解，构造了一个近似解，unCi i
i 1
利用J (Un)
11
a(Un,Un) (f ,Un)
22i
n
a(
hi
式中i,j表示节点(i,j)上的网函数。
令n(xi,yj)
ijf
『n(Xi, yj)fij
f (x, yj)
利用Taylor
展式有
(x1, yj)
2 (x,yj)
(x1,yj)
2(x, yj)
h2
4(Xi,yj)
h；
6(x, yj)
6\
h2
2
X
12
4
X
360
6
X
(h )
(xy J
2 (xy)
(Xi,yj 1)
i 1i 1
Galerkin法为
n
Rit2-Galerkin法方程：a(i,j)^( f,j)
i 1
5.有限元法：将偏微分方程转化为变分形式，选定单元的形状，对求解域作剖分，进而构造基函数或单元形状函数，形成有限元空间，将偏微分方程转化成了有限元方程，利用有效的有限元方程的解法，给出偏微分方程近似解的过程称为有限元法。
利用中心差商公式
1
[(—)](i,j)
i 1,j
i2)i，j
2
i2i 1，j
2
h
[
i,j 12i,j i,j 1
h1
将②
③两式代入①式得
I1
i
即poission
.1 i 1,j( .1
iii
2 2
h
1)i,j
2
2
1,j
r^] f(
方程极坐标形式的差分方程。
14•解：将u：+1,u：-1,u：+1, u：,按照Taylor在u：处展开整理得到其截断误差为
1.课本p2有证明
2.课本P8,P2有说明
3.课本P15, P20有说明
4.Rit2法，设Un是U的n维子空间，1,2...n是Un的一组基底，Un中的任一元素Un可
n11nn
表为UnCi i，则J(Un) -a(Un,Un)(仁叫)-a(i，jgC」Cj(f ,j)是
i 122i, j 1j 1
C1,的二次函数，a(i,j) a(j,i)，令J(Un)0，从而得到C| ,满足Cj
24t4
(Xj,
1tn)—
①
3
144
24t4
3
t3
(Xj,
2tn)
②
3
2
t3
（省去了
2
的商阶无穷小）
3
2
t
12
ah
12
44
a
12h
从而得到…(h)「(2h2)订)
15.用Fourier方法讨论向前差分格式的稳定性。
解：向前差分格式u：+=ru：++(1- 2r)u：+ru：i。以u：=vkexp(ia jh)代 入得
aj
)(N
1,juj 1
hq(x
hq
hjq(Xj 1
1
hj10逖
U1,U2,…，Un的线性方程组!
)Ui
bj
)U
)(1
)]d
h
h
h
)]d
11
0[hj 1p(Xj
h
2
h
)d
7.矩形剖分假定区域C1可以分割成为有限个互不重叠的矩形的和, 标轴平行，任意两个矩形或者不相交或者有公共的边和公共的顶点，
且每个小矩形的边和坐
,j 1
n
i,j)cCjCj( f,
j 1
j)确定G,求得近似解
Un的
过程
Galerki n
n
法：为求得UnCi i
i 1
形式
的近似解，在系数
Ci使Un关于VUn，
满足
a(Un,V)
(f,V)，对任
意
VUn或(
取Vj,1 j n
)
nn
a(i,j)Ci(f,j), j 1,2...n的情况下确定c，从而得到近似解unCi i的过程称
6.
之间的小区间Ii［Xii,xJ,hiXi人i，由节点上的一组值Uo0,Ui,U2...U|，按线
性插值公式
h
X经Ui0X Ii,i hi
1,2.确定试探空间
Un，令
F(x)
把1i变到
轴上的参考但愿［0,1］令No()1
,N
则:
U
N0( )ui 1
N
)Ui,x Ii0将®带入该函数
J(u)
1
2
2
qu
分割成三角形网，称为G的三角剖分。
基函数的取法：通过构造Lagrange型插值公式可以得到基函数的取法。不妨以R(x,y)是
一次多项式为例，得到R(x, y) L!1L22L33，其中U是相应于节点1的基函数在 △上的限制(具体的过程，可参考课本：P57P58)
9.题，参考课后习题P92的第一题，具体过程可参考积分插值的推导过程
2fu)dx
得到:
J(U
b
a
2
qu
2 fu )dx
2i
I
qu
h fUndX
带入C2可得
J(U
J (Un)
令uj
其中
aj
1
0
1
hi
i0
1
1,juj 1
f (X
h
ajjUj
1
0
1
0
1
0
1
b
aj 1,j
aj 1,j
aj,j
从而得到
hj
1P(X1
1
h
hj 1
)]
h
)_(U
\2
U
hi
hq(X
hi
)(N
2
)U
hi