计算地球物理学Word文档下载推荐.docx
- 文档编号:19561803
- 上传时间:2023-01-07
- 格式:DOCX
- 页数:19
- 大小:605.24KB
计算地球物理学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)
计算结果:
近似解图示
结果:
近似解图:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算 地球物理学
![提示](https://static.bdocx.com/images/bang_tan.gif)