计算地球物理学

计算地球物理学
成果报告
QQ 593223044
计算地球物理学
一、双曲线形方程算法和程序：
在如下问题中，对下列给定定值，用程序求解波动方程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<x<1,0<t<0.1,初始条件为：
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<x<4
u(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.5
u(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;
%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');
end
问题1：
（a）
P1=54.2857 p2=41.4286 p3=36.4286 p4=75.7143 p5=65.0000 p6=54.2857 p7=93.5714 p8=88.5714 p9=75.7143 （b）
计算结果：
近似解图示
问题2：结果：
近似解图：。