《地球物理计算方法重难点题库》由会员分享,可在线阅读,更多相关《地球物理计算方法重难点题库(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