《地球物理计算方法重难点题库》由会员分享，可在线阅读，更多相关《地球物理计算方法重难点题库（14页珍藏版）》请在金锄头文库上搜索。
1、地球物理计算方法重难点题库一、双曲线形方程算法和程序：在如下问题中，对下列给定定值，用程序求解波动方程ux（x，t）=c2uxx(x,t)，其中0xa且0tb，边界条件为：u(0,t)=0且u(a,t)=0, 0tbu(x,0)=f(x), 0xaux(x,0)=g(x), 0xa用surf和contour命令画图得到近似值解。1. 设a=1,b=1,c=1,f(x)=sin(x)，g(x)=0。为了方便起见，选择h=0.1,k=0.1。2. 设a=1,b=1,c=2,f(x)=x-x2，g(x)=0。为了方便起见，选择h=0.1,k=0.05。解：程序代码：function u,r ,x,y
2、 = finedif( f,g,a,b,c,h,k )%finedion 波动方程的差分方法程序% f：初始条件方程，字符型（sring）；% g：边界条件方程，字符型（sring）；% a：位置x的上限0,a；% b：时间t的上限0,b；% c：方程系数；% h：x的剖分步长；% k：t的剖分步长；n=a/h+1;m=b/k+1;r=c*k/h;r2=r2;r22=r2/2;s1=1-r2;s2=2-2*r2;U=zeros(n,m);%赋值边界条件for i=2:n-1 U(i,1)=feval(f,h*(i-1); U(i,2)=s1*feval(f,h*(i-1)+k*feval(g,
3、h*(i-1)+r22*(feval(f,h*i)+feval(f,h*(i-2);end%求取个点数值for j=3:m for i=2:(n-1) U(i,j)=s2*U(i,j-1)+r2*(U(i-1,j-1)+U(i+1,j-1)-U(i,j-2); end endu=U%坐标量展示：x=0:h:a;y=0:k:bend问题1：稳定性条件分析与运算结果：r=1 结果稳定结果图展示： 问题2：稳定性条件分析与运算结果：r=1 结果稳定 结果图展示：二、抛物型方程的算法和程序： 求解热传导方程：ut(x,t)=c2uxx(x,t),其中0x1,0t0.1,初始条件为：u(x,0)=f(x
4、),边界条件为：u(0,t)=g1(t),u(1,t)=g2(x)。对给定的值使用surf和contour命令画近似解。1、 使用f（x）=sin(x)+sin(2x), g1(x)=g2(x)=0, h=0.1,k=0.005.2、 使用f（x）=3-|3x-1|-|3x-2|, g1(x)=t2 ,g2(x)=e, h=0.1,k=0.005.解：程序代码：function u,r,x,y = forwdif(f,g1,g2,a,b,c,h,k )% forwdif抛物线型方程的解法% f：初始条件，字符型（string）；% g1，g2：左右边界条件，字符型（string）；% a：位置
5、上限0,a；% b：时间上线0,b；% c：方程系数；% h,k：位置和时间的剖分步长； n=a/h+1;m=b/k+1;r=c2*k/h2;s=1-2*r;U=zeros(n,m);%赋值边界条件 U(n,1:m)=feval(g2,0:k:b);U(1,1:m)=feval(g1,0:k:b);%赋值初始条件U(2:n-1,1)=feval(f,h:h:(n-2)*h);%计算for j=2:m for i=2:n-1 U(i,j)=s*U(i,j-1)+r*(U(i-1,j-1)+U(i+1,j-1); end endendu=U问题1：稳定性条件分析与运算结果：r=0.5 结果稳定结果
6、图展示：问题2：稳定性条件分析与运算结果：r=0.5 结果稳定 三、椭圆形方程算法和程序：1、（a）用程序计算5*5的网格，确定9个未知数平p1、p2p9的方程组，来求解矩形区域R=(x,y)|0x4,0y4内的谐波函数u(x,y)的近似值。,边界为：u(x,0)=10和u(x,4)=120,0x4u(0,y)=90和u(4,y)= 40,0y4（b）用9*9的网格求解近似解。2、用程序计算矩形区域R=(x,y)|0x1.5,0y1.5内的谐波函数u(x,y)的近似值，h=0.15,边界为：u(x,0)=x4和u(x,4)=x4-13.5x2+5.0625,0x1.5u(0,y)=y4和u(4
7、,y)=y4-13.5y4+5.0625,0ytol) & (cnt=max1) err=0; for j=2:m-1 for i=2:n-1 relx=w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+U(i-1,j)-4*U(i,j)/4; U(i,j)=U(i,j)+relx; if (err=abs(relx) err=abs(relx); end end endcnt=cnt+1;endu=flipud(U);end问题1：（a）P1=54.2857 p2=41.4286 p3=36.4286p4=75.7143 p5=65.0000 p6=54.2857 p7=93.5714 p8=88.5714 p9=75.7143（b）计算结果：近似解图示问题2：结果： 近似解图：14
