各向同性介质地震波场模拟，速度-应力弹性波方程交错网格有限差分法

基础源码如下，基本一看就懂，没有复杂的东西，也没有边界吸收啥的，非常适合初学者快速入门（希望读者阅读后都能觉得原来波动方程数值模拟竟如此简单，哈哈）。进阶版（考虑各向异性，PML边界等）后面有更新的帖子，欢迎关注收藏，不懂的随时交流。

#include<stdio.h>   
#include<math.h>   
#include<iostream>   
#include<fstream>   
#include<iomanip>   
using namespace std;

#define PI 3.1415926   
#define dt 0.0005   
#define dx 10   
#define dz 10   
#define NX 401   
#define NZ 401   
#define NT 1201   
#define N 1 

float Txx[NX][NZ], Tzz[NX][NZ], Txz[NX - 1][NZ - 1], Vx[NX - 1][NZ], Vz[NX][NZ - 1], L[NX][NZ], M[NX][NZ], e[NX][NZ];

int main()
{
	int i,j,k,f0=10;float t0=1.2/f0;
	FILE *fp1,*fp2,*fp3,*fp4,*fp5;
	for(i=0;i<NX;i++)
		for(j=0;j<NZ;j++)
		{
			Txx[i][j]=0;
			Tzz[i][j]=0;
			//L[i][j]=1.19*pow(10, 10);
			L[i][j]=0.19*pow(10, 10);
			M[i][j]=5.4*pow(10, 9);
			e[i][j]=2.0*pow(10, 3);
		}
	for(i=0;i<NX-1;i++)
		for(j=0;j<NZ;j++)
		{
			Vx[i][j]=0;
		}
	for(i=0;i<NX;i++)
		for(j=0;j<NZ-1;j++)
		{
			Vz[i][j]=0;
		}
	for(i=0;i<NX-1;i++)
		for(j=0;j<NZ-1;j++)
		{
			Txz[i][j]=0;
		}


	for(k=0;k<NT;k++)
	{
		for(i=N;i<NX-N;i++)
			for(j=N;j<NZ-N;j++)
			{
				Txx[i][j]=Txx[i][j]+dt*((L[i][j]+2*M[i][j])*(Vx[i][j]-Vx[i-1][j])/dx+(Vz[i][j]-Vz[i][j-1])/dz*L[i][j]);
				Tzz[i][j]=Tzz[i][j]+dt*((L[i][j]+2*M[i][j])*(Vz[i][j]-Vz[i][j-1])/dz+(Vx[i][j]-Vx[i-1][j])/dx*L[i][j]);
				if(i==NX/2&&j==NZ/2)//这里加雷克子波震源
				{
					Txx[i][j]=Txx[i][j]+1000*(1-2*PI*f0*(k*dt-t0)*PI*f0*(k*dt-t0))*exp(-PI*f0*(k*dt-t0)*PI*f0*(k*dt-t0));
					Tzz[i][j]=Tzz[i][j]+1000*(1-2*PI*f0*(k*dt-t0)*PI*f0*(k*dt-t0))*exp(-PI*f0*(k*dt-t0)*PI*f0*(k*dt-t0));
				}
			}

		for(i=N;i<NX-N-1;i++)
			for(j=N;j<NZ-N-1;j++)
				Txz[i][j]=Txz[i][j]+dt*M[i][j]*((Vz[i+1][j]-Vz[i][j])/dx+(Vx[i][j+1]-Vx[i][j])/dz);

		for(i=N;i<NX-N-1;i++)
			for(j=N;j<NZ-N;j++)
			{
				Vx[i][j]=Vx[i][j]+dt/e[i][j]*((Txx[i+1][j]-Txx[i][j])/dx+(Txz[i][j]-Txz[i][j-1])/dz);
				//if(i==NX/2&&j==NZ/2)
				//{
				//	Vx[i][j]=Vx[i][j]+1000*(1-2*PI*f0*(k*dt-t0)*PI*f0*(k*dt-t0))*exp(-PI*f0*(k*dt-t0)*PI*f0*(k*dt-t0));
				//}
			}

		for(i=N;i<NX-N;i++)
			for(j=N;j<NZ-N-1;j++)
			{
				Vz[i][j]=Vz[i][j]+dt/e[i][j]*((Tzz[i][j+1]-Tzz[i][j])/dz+(Txz[i][j]-Txz[i-1][j])/dx);
				//if(i==NX/2&&j==NZ/2)
				//{
				//	Vz[i][j]=Vz[i][j]+1000*(1-2*PI*f0*(k*dt-t0)*PI*f0*(k*dt-t0))*exp(-PI*f0*(k*dt-t0)*PI*f0*(k*dt-t0));
				//}
			}

	


	}
		fp1=fopen("Txx.dat", "wb+");
		for (i = 0; i < NX; i++)
			for (j = 0; j < NZ; j++)
				fwrite(&Txx[i][j], 4, 1, fp1);
		fp2=fopen("Tzz.dat", "wb+");
		for (i = 0; i < NX; i++)
			for (j = 0; j < NZ; j++)
				fwrite(&Tzz[i][j], 4, 1, fp2);
		fp3=fopen("Txz.dat", "wb+");
		for (i = 0; i < NX-1; i++)
			for (j = 0; j < NZ-1; j++)
				fwrite(&Txz[i][j], 4, 1, fp3);
		fp4=fopen("Vx.dat", "wb+");
		for (i = 0; i < NX-1; i++)
			for (j = 0; j < NZ; j++)
				fwrite(&Vx[i][j], 4, 1, fp4);
		fp5=fopen("Vz.dat", "wb+");
		for (i = 0; i < NX; i++)
			for (j = 0; j < NZ-1; j++)
				fwrite(&Vz[i][j], 4, 1, fp5);
}

下图中的FSM和SPAFSM指的是走时场的分布，是波动方程的高频近似：程函方程计算得到的，可以参考快速扫描法FSM方法的帖子。