计算地球物理学.docx
- 文档编号:6548044
- 上传时间:2023-01-07
- 格式:DOCX
- 页数:19
- 大小:605.24KB
计算地球物理学.docx
《计算地球物理学.docx》由会员分享,可在线阅读,更多相关《计算地球物理学.docx(19页珍藏版)》请在冰豆网上搜索。
计算地球物理学
计算地球物理学
成果报告
QQ593223044
计算地球物理学
一、双曲线形方程算法和程序:
在如下问题中,对下列给定定值,用程序求解波动方程ux(x,t)=c2uxx(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
ux(x,0)=g(x),0≤x≤a
用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);
%赋值边界条件
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
end
u=U'
%坐标量展示:
x=0:
h:
a;
y=0:
k:
b
end
问题1:
稳定性条件分析与运算结果:
r=1结果稳定
结果图展示:
问题2:
稳定性条件分析与运算结果:
r=1结果稳定
结果图展示:
二、抛物型方程的算法和程序:
求解热传导方程:
ut(x,t)=c2uxx(x,t),其中0 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抛物线型方程的解法¨ %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)'; %计算 forj=2: m fori=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; %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 fori=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.2857p2=41.4286p3=36.4286 p4=75.7143p5=65.0000p6=54.2857 p7=93.5714p8=88.5714p9=75.7143 (b) 计算结果: 近似解图示 问题2: 结果: 近似解图:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算 地球物理学