计算地球物理学

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

一、双曲线形方程算法和程序：

在如下问题中，对下列给定定值，用程序求解波动方程u x（x，t）=c2u xx(x,t)，其中0≤x≤a且0≤t≤b，边界条件为：

u(0,t)=0且u(a,t)=0, 0≤t≤b

u(x,0)=f(x), 0≤x≤a

(x,0)=g(x), 0≤x≤a

u

x

用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] = 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);

%赋值边界条件

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,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

end

u=U'

%坐标量展示：

x=0:h:a;

y=0:k:b

end

问题1：

稳定性条件分析与运算结果：

r=1 结果稳定

结果图展示：

问题2：

稳定性条件分析与运算结果：r=1 结果稳定

结果图展示：

二、抛物型方程的算法和程序：

求解热传导方程：u

t (x,t)=c2u

xx

(x,t),其中0

u(x,0)=f(x),边界条件为：u(0,t)=g

1(t),u(1,t)=g

2

(x)。对给定的值使用surf

和contour命令画近似解。

1、使用f（x）=sin(πx)+sin(2πx), g

1(x)=g

2

(x)=0, h=0.1,k=0.005.

2、使用f（x）=3-|3x-1|-|3x-2|, g

1(x)=t2 ,g

2

(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)';

%计算

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

end

end

u=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

u(0,y)=90和u(4,y)= 40,0

（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

u(0,y)=y4和u(4,y)=y4-13.5y4+5.0625,0

程序代码：

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;

%SOR parameter（超松弛因子）

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;

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

end

cnt=cnt+1;

end

u=flipud(U');