计算地球物理学Word文档下载推荐.docx
《计算地球物理学Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《计算地球物理学Word文档下载推荐.docx（19页珍藏版）》请在冰豆网上搜索。
程序代码：
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-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
%求取个点数值
forj=3:
m
fori=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）,其中0<
x<
1,0<
t<
0.1,初始条件为：
u（x,0）=f（x）,边界条件为：
u（0,t）=g1（t）,u（1,t）=g2（x）。
对给定的值使用surf和contour命令画近似解。
1、使用f（x）=sin（πx）+sin（2πx）,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：
位置和时间的剖分步长；
r=c^2*k/h^2;
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）'
;
%计算
forj=2:
U（i,j）=s*U（i,j-1）+r*（U（i-1,j-1）+U（i+1,j-1））;
end
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<
4
u（0,y）=90和u（4,y）=40,0<
y<
（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<
1.5
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）;
%边界条件赋值
m）=feval（f3,0:
（m-1）*h）'
m）=feval（f4,0:
U（1:
n,1）=feval（f1,0:
（n-1）*h）;
n,m）=feval（f2,0:
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-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.2857p2=41.4286p3=36.4286
p4=75.7143p5=65.0000p6=54.2857
p7=93.5714p8=88.5714p9=75.7143
（b）
计算结果：
近似解图示
结果：
近似解图：
升级会员 