计算地球物理学(2)

二、抛物型方程的算法和程序：

求解热传导方程：ut(x,t)=c2uxx(x,t),其中0

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)'; %计算

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 结果稳定