计算地球物理学Word文档下载推荐.docx

1、程序代码：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=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*feva

2、l（f,h*（i-1）+k*feval（g,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 u=U%坐标量展示：x=0:h:a;y=0:k:b问题1：稳定性条件分析与运算结果：r=1 结果稳定结果图展示：问题2：二、抛物型方程的算法和程序： 求解热传导方程：ut（x,t）=c2uxx（x,t）,其中0x1,0t0.1,初始条件为：u（x,0）=f（x）,边界条件为：u（0,t

3、）=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抛物线型方程的解法初始条件，字符型（string）；% g1，g2：左右边界条件，字符型（string）；位置上限0,a；时间上线0,b；% h,k：位置和时间的剖分

4、步长；r=c2*k/h2;s=1-2*r;%赋值边界条件 U（n,1:m）=feval（g2,0:b）;U（1,1:m）=feval（g1,0:赋值初始条件U（2:n-1,1）=feval（f,h:（n-2）*h）;%计算for j=2: U（i,j）=s*U（i,j-1）+r*（U（i-1,j-1）+U（i+1,j-1）; endr=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,04u（0,y）

5、=90和u（4,y）= 40,0y（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,01.5u（0,y）=y4和u（4,y）=y4-13.5y4+5.0625,0tol） & （cnt=max1） err=0; for j=2:m-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）;cnt=cnt+1;u=flipud（U）;（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）计算结果：近似解图示结果：近似解图：