计算地球物理学
成果报告
QQ 593223044
www.zggxbbs.com
计算地球物理学
一、双曲线形方程算法和程序：
在如下问题中，对下列给定定值，用程序求解波动方程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); %赋值边界条件
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 结果稳定
结果图展示：