地球物理计算方法重难点题库

地球物理计算方法重难点题库一、双曲线形方程算法和程序：在如下问题中，对下列给定定值，用程序求解波动方程ux（x，t）=c2uxx(x,t)，其中0≤x≤a且0≤t≤b，边界条件为：u(0,t)=0且u(a,t)=0,0≤t≤bu(x,0)=f(x),0≤x≤aux(x,0)=g(x),0≤x≤a用surf和contour命令画图得到近似值解。设a=1,b=1,c=1,f(x)=sin(πx)，g(x)=0。为了方便起见，选择h=0.1,k=0.1。设a=1,b=1,c=2,f(x)=x-x2，g(x)=0。为了方便起见，选择h=0.1,k=0.05。解：程序代码：function[u,r,x,y]=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=r^2;r22=r^2/2;s1=1-r^2;s2=2-2*r^2;U=zeros(n,m);%赋值边界条件fori=2:n-1U(i,1)=feval(f,h*(i-1));U(i,2)=s1*feval(f,h*(i-1))+k*feval(g,h*(i-1))+r22*(feval(f,h*i)+feval(f,h*(i-2)));end%求取个点数值forj=3:mfori=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);endendu=U'%坐标量展示：x=0:h:a;y=0:k:bend问题1：稳定性条件分析与运算结果：r=1结果稳定结果图展示：问题2：稳定性条件分析与运算结果：r=1结果稳定结果图展示：二、抛物型方程的算法和程序：求解热传导方程：ut(x,t)=c2uxx(x,t),其中0<x<1,0<t<0.1,初始条件为：u(x,0)=f(x),边界条件为：u(0,t)=g1(t),u(1,t)=g2(x)。对给定的值使用surf和contour命令画近似解。使用f（x）=sin(πx)+sin(2πx),g1(x)=g2(x)=0,h=0.1,k=0.005.使用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：位置上限[0,a]；%b：时间上线[0,b]；%c：方程系数；%h,k：位置和时间的剖分步长；n=a/h+1;m=b/k+1;r=c^2*k/h^2;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)';%计算forj=2:mfori=2:n-1U(i,j)=s*U(i,j-1)+r*(U(i-1,j-1)+U(i+1,j-1));endendendu=U'问题1：稳定性条件分析与运算结果：r=0.5结果稳定结果图展示：问题2：稳定性条件分析与运算结果：r=0.5结果稳定三、椭圆形方程算法和程序：1、（a）用程序计算5*5的网格，确定9个未知数平p1、p2……p9的方程组，来求解矩形区域R={(x,y)|0≤x≤4,0≤y≤4}内的谐波函数u(x,y)的近似值。,边界为：u(x,0)=10和u(x,4)=120,0<x<4u(0,y)=90和u(4,y)=40,0<y<4（b）用9*9的网格求解近似解。2、用程序计算矩形区域R={(x,y)|0≤x≤1.5,0≤y≤1.5}内的谐波函数u(x,y)的近似值，h=0.15,边界为：u(x,0)=x4和u(x,4)=x4-13.5x2+5.0625,0<x<1.5u(0,y)=y4和u(4,y)=y4-13.5y4+5.0625,0<y<1.5程序代码：function[u,w,x,y]=dirich(f1,f2,f3,f4,a,b,h,tol,max1)%function：椭圆型方程的差分法%f1,f2,f3,f4：边界条件，字符型（string）%a，b：x，y上限值；%h：步长n=fix(a/h)+1;m=fix(b/h)+1;ave=(a*(feval(f1,0)+feval(f2,0))+b*(feval(f3,0)+feval(f4,0)))/(2*a+2*b);U=ave*ones(n,m);%边界条件赋值U(1,1:m)=feval(f3,0:h:(m-1)*h)';U(n,1:m)=feval(f4,0:h:(m-1)*h)';U(1:n,1)=feval(f1,0:h:(n-1)*h);U(1:n,m)=feval(f2,0:h:(n-1)*h);U(1,1)=(U(1,2)+U(2,1))/2;U(1,m)=(U(1,m-1)+U(2,m))/2;U(n,1)=(U(n-1,1)+U(n,2))/2;U(n,m)=(U(n-1,m)+U(n,m-1))/2;%SORparameter（超松弛因子）w=4/(2+sqrt(4-(cos(pi/(n-1))+cos(pi/(m-1)))^2));%计算err=1;cnt=0;while((err>tol)&(cnt<=max1))err=0;forj=2:m-1fori=2:n-1relx=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);endendendcnt=cnt+1;endu=flipu