下载文档原格式
// ParallelPI.cpp : 定义控制台应用程序的入口点。
//// ParallelPI.cpp : 定义控制台应用程序的入口点。
////地层为三层速度模型，使用的进程数为4个，零进程为主进程，负责接收计算结果，其他进程进行计算#include "mpi.h"#include <stdio.h>#include <stdlib.h>#include <math.h>void main(int argc, char* argv[]){FILE *fp;double pi, r, f;int sourcex;int sourcey;int MAX_Column;//列数int MAX_Row;//矩阵的行数int myid, numprocs, master, column, row, sourceid, T,sendnum,recvnum;int i, j, k, n, position;//计数变量int up, down;int begin, end;//定义开始结束计算的行数double dx, dz, dt;//时空间步长double ***u;//时间跨度三层double **v;//速度矩阵double s[400], **temp;//扰动MPI_Status status;MPI_Init(&argc, &argv);MPI_Comm_size(MPI_COMM_WORLD, &numprocs);//numprocs=4MPI_Comm_rank(MPI_COMM_WORLD, &myid);//给定正演模拟的参数pi = 3.1415926, r = 4.0, f = 25.0;sourcex = 150;//震源位置sourcey = 150;MAX_Column = 300;MAX_Row = 300;column = MAX_Column / (numprocs - 1);row = MAX_Row;begin = 2, end = column + 1;T = 400;dx = 5.0, dz = 5.0, dt = 0.001;for (k = 0; k < T; k++) { s[k] = exp(-pow((2 * pi*f*k*dt / r), 2))*cos(2 * pi*f*k*dt); }//定义扰动，雷克子波//定义速度矩阵v = (double **)malloc(sizeof(double *)*(column + 4));for (i = 0; i < column + 4; i++){v[i] = (double *)malloc(sizeof(double)*row);}if (myid == 0){if ((fp = fopen("C:\\Users\\User\\Desktop\\velocity.txt", "r")) != NULL){for (sendnum = 1; sendnum <= 3; sendnum++){for (i = 2; i <= column + 1; i++){for (j = 0; j < row; j++) {fscanf(fp, "%lf", &v[i][j]);}}for (i = 2; i <= column + 1; i++)MPI_Send(v[i], row, MPI_DOUBLE, sendnum, 1, MPI_COMM_WORLD);}fclose(fp);}}if (myid != 0){for ( i=2; i<=column+1; i++)MPI_Recv(v[i], row, MPI_DOUBLE, 0, 1, MPI_COMM_WORLD, &status);}//指定每一个进程开始和结束计算的行号if (myid == 1){for (j = 0; j<row; j++)//v[2][j]=0.0;begin = 4;}if (myid == numprocs - 1){for (j = 0; j<row; j++)//v[column+1][j]=0.0;end = column - 1;}//u记录各点振幅值u = (double ***)malloc(sizeof(double **)* 3);for (k = 0; k<3; k++){u[k] = (double **)malloc(sizeof(double *)*(column + 4));for (i = 0; i<column + 4; i++){u[k][i] = (double *)malloc(row*sizeof(double));}}for (k = 0; k<3; k++)//赋初值for (i = 0; i<column + 4; i++)for (j = 0; j<row; j++)u[k][i][j] = 0;if (myid>1)//定义左右进程虚进程up = myid - 1;elseif (myid != 0) up = MPI_PROC_NULL;if (myid<numprocs - 1 && myid != 0)down = myid + 1;elseif (myid != 0) down = MPI_PROC_NULL;sourceid = sourcey / column + 1;//确定震源应该位于的进程//进程间通信并且进行计算if (myid != 0)for (k = 0; k<T; k++){MPI_Sendrecv(u[1][column], row, MPI_DOUBLE, down, 1, u[1][0], row, MPI_DOUBLE, up, 1, MPI_COMM_WORLD, &status);//从上向下平移MPI_Sendrecv(u[1][column + 1], row, MPI_DOUBLE, down, 1, u[1][1], row, MPI_DOUBLE, up, 1, MPI_COMM_WORLD, &status);//从上向下平移MPI_Sendrecv(u[1][2], row, MPI_DOUBLE, up, 2, u[1][column + 2], row, MPI_DOUBLE, down, 2, MPI_COMM_WORLD, &status);//从下向上平移MPI_Sendrecv(u[1][3], row, MPI_DOUBLE, up, 2, u[1][column + 3], row, MPI_DOUBLE, down, 2, MPI_COMM_WORLD, &status);//从下向上平移for (i = begin; i <= end; i++){for (j = 2; j<row - 2; j++){u[2][i][j] = pow(v[i][j] * dt / dx, 2)*(-1.0 / 12 * (u[1][i + 2][j] + u[1][i - 2][j]) + 4.0 / 3 * (u[1][i + 1][j] + u[1][i - 1][j]) - 5.0 / 2 * u[1][i][j]) + pow(v[i][j] * dt / dz, 2)*(-1.0 / 12 * (u[1][i][j + 2] + u[1][i][j - 2]) + 4.0 / 3 * (u[1][i][j + 1] + u[1][i][j - 1]) - 5.0 / 2 * u[1][i][j]) + 2 * u[1][i][j] - u[0][i][j];//时间二阶精度差分格式if (myid == sourceid&&i == (sourcey - column*(myid - 1) + 2) && j == sourcex)//如果某个进程计算到的网格点正好是震源点，则应在该点的振幅值上加上原始扰动的振幅值u[2][i][j] += s[k];}}temp = u[0];u[0] = u[1];u[1] = u[2];u[2] = temp;}if (myid != 0){for (i = 2; i <= column + 1; i++)MPI_Send(u[1][i], row, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);//按行向主进程发送运算结果}//printf("%d ok", myid);//输出数据，前两行为行数及列数，后面为波前快照的数据体。
声波方程数值模拟实验报告.基础理论知识需要的已知条件包括：1）震源函数地层速度（波速） 边界条件.2.2.2苇=v 2(諾二2)S(t)一 t :x :zv （x, Z ）是介质在点（x , z ）处的纵波速度，u 为描述速度位或者压力的波场，s （t ）为震 源函数。
为求式（4-1）的数值解，必须将此式 离散化，即用有限差分来逼近导数，用差商代替 微商。
为此，先把空间模型网格化（如图4-1所示）。
设x 、z 方向的网格间隔长度为h ，厶t 为时间采样步长，则有:X 二i^h（i 为正整数） z = j.）h（j 为正整数）^n ：t（n 为正整数）u ：j 表示在（i,j ）点，k 时刻的波场值。
k1. 2)3)>2u 2f cu〒=V p ^~2 t :x22w 2 w 厂二 V （—T :t ；x 声波方程的有限差分法数值模拟 对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述:2. 弹性波方程:二b s (t) 一 z22_w )；z 2)(4-1)将u i,j在（i,j）点k时刻用Taylor展式展开：k 1 kjUu i,j 7,j ■—ct将U i k j 」在(i,j)点k 时刻用Taylor 展式展开:k k4 k k5 k[U i2jui 2,j] 3[ui 4,j ui ・1,j] —?U i,j }S(t)*、(i -i °)**「(j - j °)(4-7)式中v(i, j)为介质速度的空间离散值，：h 是空间离散步长，=t 为时间离散步长，s(k)为震源函数，关于 s(k) 一般使用一个理论的雷克型子波代替，即:上式中，t 为时间，f 为中心频率，一般取为20-40HZ ， 为控制频带宽度的参数,(4-2)k 1k ；'UUi,jUi,2*「7 2 ；:t 2(4-3)将上两式相加，略去高阶小量, 整理得(i,j)点k 时刻的二阶时间微商为:2k 1kk J；：u u i,j -2u i,j u i,j.:t 242(4-4)对于空间微分，采用四阶精度差分格式,(以X 方向为例)即将U i*；j 、*Tj 分别在(i,j)点k 时刻展开到四阶小量，消除四阶小量并解出二阶微分得:-2~ u 11 k k 4 k k 5 k{—T7[u i 二 j +u id2,j ]+：[u i4j +U i*,j ] —=u i,j } 12 3 2:x L X 2(4-5)同理可得:2；=u 1.1 kk4 kk5 k一 2 人 2{一 石[Ui,j ，+Ui,j~2]+：[Ui,j 」+Ui,j^]—；U i,j}:z-z12 32(4-6)这就实现了用网个点波场值的差商代替了偏微分方程的微商，将上三个式子代入 (4-1)式中得:k 1 Ui,j kk 4= 2u i,j —Ui,jV :-1 1 k k 4 k k—{p [u=j U i/ yij —j：h 25 k H-U i,j }h 2 {12s (t )二 e(-2 f / )2t2cos2 二 ft(4-8)般取3-5。
超声波声场的数值模拟及分析研究I. 前言超声波技术作为一种重要的非破坏性检测手段，广泛应用于医学、工业、军事等领域。
在这些领域中，对于超声波传播和反射规律的研究，对提高超声波技术检测的灵敏度和准确度至关重要。
声场是一个重要的研究方向。
声场的数值模拟和分析可以帮助研究者更好地了解超声波在不同介质中的传播规律。
II. 超声波声场模拟的基本原理(一) 声场模拟的定义声场模拟是指基于声波理论对声场进行数值计算和分析，以掌握声场分布的规律性，实现声波信号的处理和应用。
(二) 超声波传播的数学模型超声波传播的数学模型是典型的波动方程，其中声压 $ P(x,t) $ 是变量，表示声波在空间位置 $x$ 和时间 $t$ 处的强度。
根据波动方程，超声波传播的速度取决于介质的密度和压力。
(三) 超声波声场模拟方法声场模拟的方法主要有有限差分法、边界元法和声束追踪法。
其中，有限差分法适用于平面和轴对称的声场计算，边界元法适用于无限空间中的声场计算，声束追踪法适用于计算弥散源下的声场传播。
III. 超声波声场模拟的应用及优点(一)超声波生物医学应用中的声场模拟超声波在生物医学应用中广泛应用，如超声心动图、超声诊断、超声治疗等。
声场模拟可以辅助医生分析生物组织中声波的传播规律，提高检测准确度和效率，并为医疗维护提供重要支持。
(二) 超声波工业检测中的声场模拟超声波检测在工业生产中应用广泛，如对金属、塑料等材料的瑕疵检测。
声场模拟技术可以模拟材料内部的声场分布，定位和评估瑕疵，提高测试的准确性和可靠性。
(三) 超声波声场模拟的优点超声波声场模拟可以模拟不同介质中的声波传播规律，为声波技术的研究提供了可靠依据。
与实验方法相比，声场模拟不受时间和空间限制，更加灵活方便。
声场模拟还可以降低实验成本，减少实验过程中的危险和损失。
IV. 超声波声场模拟的发展趋势(一)模拟软件的发展超声波声场模拟需要强大的计算机算力和复杂的程序设计技术，目前市场上已经有多种声场模拟软件，如ABAQUS、COMSOL、ANSYS等。
一、实验目的1. 理解声波的基本特性及其传播规律。
2. 掌握使用驻波法测量声速的方法。
3. 学习使用示波器观察和分析声波信号。
4. 通过实验，提高动手能力和数据处理能力。
二、实验原理声波是一种机械波，它通过介质中的质点振动传播。
声速是指声波在介质中传播的速度，其大小与介质的性质有关。
本实验采用驻波法测量声速，其原理如下：当声波在两束反射波相遇时，会产生干涉现象，形成驻波。
驻波的特点是波节和波腹的位置固定不变，波节处的振幅为零，波腹处的振幅最大。
根据驻波波节和波腹的位置关系，可以计算出声波的波长，从而求得声速。
三、实验仪器与器材1. 实验仪器：驻波实验装置、示波器、信号发生器、频率计、测量尺等。
2. 实验器材：玻璃管、橡皮膜、连接线、电阻等。
四、实验步骤1. 将玻璃管水平放置，将橡皮膜固定在管口。
2. 将信号发生器的输出信号接入驻波实验装置，调节频率使驻波形成。
3. 打开示波器，将示波器的输入通道连接到驻波实验装置的输出端。
4. 观察示波器屏幕上的波形，找到波节和波腹的位置。
5. 使用测量尺测量波节和波腹之间的距离，得到声波的波长。
6. 使用频率计测量信号发生器的输出频率，得到声波的频率。
7. 根据声波的波长和频率，计算声速。
五、实验数据1. 波长测量数据：波节位置：L1 = 5.5 cm波腹位置：L2 = 15.6 cm波长：λ = L2 - L1 = 10.1 cm = 0.101 m2. 频率测量数据：频率：f = 500 Hz3. 声速计算数据：声速：v = λf = 0.101 m × 500 Hz = 50.5 m/s六、实验结果与分析通过实验，我们测量得到了声波的波长和频率，进而计算出了声速。
实验结果显示，在实验条件下，声速约为50.5 m/s。
与理论值相比，实验结果存在一定的误差，这可能是由于实验器材的精度限制、实验环境的影响以及实验操作中的误差等因素造成的。
七、实验结论1. 通过本实验，我们掌握了使用驻波法测量声速的方法。
声波方程有限差分数值模拟的变网格步长算法声波方程是描述声波在介质中传播的方程，在进行数值模拟时需要使用有限差分方法来近似求解。
有限差分方法将连续的空间和时间离散化，通过在离散的网格节点上计算声压场的数值来模拟声波的传播过程。
在有限差分数值模拟中，选择合适的网格步长对数值结果的精度和计算效率具有重要影响。
变网格步长算法是一种能够根据需要在不同区域自适应地调整网格步长的技术。
在声波方程有限差分数值模拟中，声波会在不同介质中以不同的速度传播，因此在网格步长选择上需要考虑介质的变化。
当声波在介质变化剧烈的区域传播时，使用较小的网格步长可以提高模拟的精度。
而在介质变化缓和的区域，使用较大的网格步长能够减少计算量。
变网格步长算法的基本思想是通过对声波传播区域进行划分，然后根据介质的变化情况调整不同区域的网格步长。
具体步骤如下：1.对传播区域进行划分：根据介质的变化情况，将声波传播区域划分为多个区域，每个区域具有不同的网格步长。
2.确定初始网格步长：在初始时，可以根据经验或者初步模拟结果选择一个合适的网格步长。
3.计算声波传播：在每个区域内分别进行声波方程的有限差分数值模拟，使用当前区域的网格步长进行计算。
4.判断误差与精度：通过计算得到的声压场数值和预期精度进行比较，如果达到要求，则结束模拟；否则，进行下一步。
5.调整网格步长：根据当前区域的模拟误差情况，调整该区域的网格步长，使其更适应声波在当前区域的传播特性。
可以根据误差大小和梯度等因素进行调整。
6.重复模拟：根据调整后的网格步长，重新进行声波方程的有限差分数值模拟，并返回步骤4通过以上步骤，可以实现声波方程有限差分数值模拟的变网格步长算法。
这种算法能够根据介质的变化情况自适应地调整网格步长，提高模拟的精度和计算效率。
然而，变网格步长算法的实现并不简单，需要合理设置划分区域和调整步长的策略，并在计算过程中不断优化和调整，才能达到较好的数值模拟效果。
声波方程有限差分数值模拟的变网格步长算法声波方程有限差分数值模拟是一种常用的声波传播模拟方法，可以在计算机上通过数值计算求解声波传播的过程。
在进行这种数值模拟时，常常需要选择合适的网格步长，以保证计算结果的准确性和计算效率。
本文将介绍一种变网格步长算法，用于优化声波方程有限差分数值模拟的计算。
声波方程可以用下面的形式表示：∂^2p/∂t^2=c^2∇^2p其中p是声场变量，t是时间，c是声速，∇^2是Laplace算子。
为了将声波方程用有限差分方法进行离散化计算，我们需要将空间和时间分别离散化。
首先，将空间离散化为网格，在每个网格点上计算声场的值。
其次，将时间离散化为离散的时间步长，通过迭代计算不同时间步长上的声场分布。
为了保证计算结果的准确性，网格步长应当满足Nyquist采样定理的要求。
即网格步长应小于声波的最小波长的一半。
根据声波方程的性质，我们可以通过声速和最高频率来估计声波的最小波长。
然后，我们可以根据最小波长来选择合适的网格步长。
然而，在实际的声波传播计算中，声场的变化往往不是均匀的。
有些区域的声场变化较大，而其他区域的声场变化较小。
如果我们在整个计算区域都采用较小的网格步长，将会造成计算资源的浪费。
因此，需要一种方法能够根据声场的变化情况来自适应地调整网格步长。
变网格步长算法就是一种能够根据声场变化情况自动调整网格步长的算法。
其基本思想是根据声场在不同网格上的变化率来决定每个网格上的网格步长。
具体的算法步骤如下：1.初始化：选择一个合适的初始网格步长。
通常可以选择根据声波的最小波长来确定。
2.计算网格步长：在每个时间步长上，对于每个网格点，计算其周围网格点上的声场变化率。
常用的方法是计算声场在三个相邻时间步长上的差分值，然后取绝对值并求平均。
根据声场变化率，调整当前网格点上的网格步长。
变化率大的网格点应该有更小的网格步长，而变化率小的网格点则可以有更大的网格步长。
3.更新声场：根据调整后的网格步长，更新所有网格点上的声场值。
声波波动方程正演模拟程序程序介绍：第一部分：加载震源，此处选用雷克子波当作震源。
编写震源程序后，我将输出的数据复制，然后我用excel做成了图片，以检验程序编写是否正确。
以下为雷克子波公式部分的程序：for(it=0;it<Nt;it++){t1=it*dt;t2=t1-t0;source[it]=(1.0-2.0*pow(PI*fm*t1,2.0))*exp(-pow(PI*fm*t1,2.0));fprintf(fp,"%.8f %.8f\n",t1,source[it]);}此处，为了成图完整，我用的是t2，而不是t1，也就是把雷克子波向右移动了一段距离，使主要部分都显示出来。
（频率采用的是30hz)从图中可以看出程序是正确的，符合理论上雷克子波的波形。
第二部分：主程序，编写声波正演模拟算子。
首先定义了各种变量，然后指定震源位置，选择权系数，给速度赋值，然后是差分算子的编写，这是主要部分，最后再进行时间转换，即把n-1时刻的值给n时刻，把n时刻的值给n+1时刻。
此处，我编写的是均匀介质声波方程规则网格的正演模拟程序，时间导数采用二阶中心差分、空间导数为2N阶差分精度，网格大小为200*200，总时间为400。
第三部分：这一部分就是记录文件。
首先记录Un文件，然后记录record文件。
模型构建与试算：1、我首先建立了一个均匀介质模型，首先利用不同时间，进行了数值模拟，得到波场快照如图所示：100ms 200ms 300ms此处，纵波速度为v=3000m/s。
模型大小为200×200，空间采样间隔为dx=dz=10m。
采用30Hz的雷克子波作为震源子波，时间采样间隔为1ms，图中可以看出，波场快照中的同相轴是圆形的，说明在均匀各向同性介质中，点源激发的波前面是一个圆，这与理论也是吻合的。
并且随着时间的增大，波前面的面积逐渐增大，说明地震波从震源中心向外传播。
声波方程数值模拟实验报告一．基础理论知识需要的已知条件包括：1．1）震源函数2）地层速度（波速） 3）边界条件2．弹性波方程：⎪⎪⎩⎪⎪⎨⎧∂∂+∂∂=∂∂+∂∂+∂∂=∂∂)()()(22222222222222z w x w v t w t S zu x uv t u s p声波方程的有限差分法数值模拟对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：)()(2222222t S zu x u v t u +∂∂+∂∂=∂∂ （4-1） (,)v x z 是介质在点（x , z ）处的纵波速度，u 为描述速度位或者压力的波场，)(t s 为震源函数。
为求式（4-1）的数值解，必须将此式离散化，即用有限差分来逼近导数，用差商代替微商。
为此，先把空间模型网格化（如图4-1所示）。
设x 、z 方向的网格间隔长度为h ∆，t ∆为时间采样步长，则有：h i x ∆= （i 为正整数）h j z ∆= (j 为正整数)t n t =∆ (n 为正整数)k j i u , 表示在(i,j)点，k 时刻的波场值。
将1,+k j i u 在(i,j)点k 时刻用Taylor 展式展开：)(*21*22*22*,1,t o t t ut tu uutk t t k t k ji k ji ∆+∆∂∂+∆∂∂+=∆=∆=+(4-2)将1,-k j i u 在(i,j)点k 时刻用Taylor 展式展开：)(*21*22*22*,1,t o t t u t tu uutk t tk t k ji k ji ∆+∆∂∂+∆∂∂-=∆=∆=-(4-3)将上两式相加，略去高阶小量，整理得(i,j)点k 时刻的二阶时间微商为：21,,1,222t u u u t u k ji k j i k j i ∆+-=∂∂-+ (4-4)对于空间微分，采用四阶精度差分格式，（以X 方向为例）即将1,2++k j i u 、1,1++k j i u 、1,1+-k j i u1,1++k j i u 分别在(i,j)点k 时刻展开到四阶小量，消除四阶小量并解出二阶微分得：}25][34][121{1,,1,1,2,2222k j i kj i k j i k j i k j i u u u u u xx u -+++-∆=∂∂+-+- (4-5) 同理可得：}25][34][121{1,1,1,2,2,222k j i kj i k j i k j i k j i u u u u u z z u -+++-∆=∂∂+-+- (4-6)这就实现了用网个点波场值的差商代替了偏微分方程的微商，将上三个式子代入(4-1)式中得：}25][34][121{2,,1,1,2,22221,,1,k j i k ji k j i k j i k j i k ji k ji k ji u u u u u h t v uuu-+++-∆∆+-=+-+--+k j i kj i k j i k j i k j i u u u u u h t v ,,1,1,2,222225][34][121{-+++-∆∆++-+-})(**)(*)(00j j i i t s --+δδ(4-7)式中),(j i v 为介质速度的空间离散值，h ∆是空间离散步长，t ∆为时间离散步长，)(k s 为震源函数，关于)(k s 一般使用一个理论的雷克型子波代替，即：ft t f t s e πγπ2cos )/2()(22-= （4-8）上式中，t 为时间， f 为中心频率，一般取为20-40HZ ，γ为控制频带宽度的参数，一般取3-5。
在实际计算过程中，需把此震源函数离散，参与波场计算。
)(*)(00j j i i --δδ确定震源位置。
稳定性条件：83/*max <∆∆h t v (4-8)这里m ax v 表示的是地下介质的最大波速；若地下介质网格间隔、最小速度、及时间采样间隔不符合(4-8)式时，第推求解（4-7）式，波场值会出现误差（高阶小量）累积，出现不稳定现象。
频散关系式：)/(min N Gf v h ≤∆ (4-9)式中m in v 为最小速度，N f 为Nyquist 频率。
一般取震源子波中的主频f 的2倍值参与计算，G 为每个波长所占的网格点数，对于空间二阶差分、时间二阶差分G 取8，而对于空间为四阶差分的情况则G 取4方能有效减少频散。
二．实验步骤1、 应用声波方程作为正演模拟的波动方程，忽略转换波的产生、传播；2、 将所提供震源函数离散后绘图；震源函数为雷克子波，离散绘图如下： fm=30; r=3; t=0.002; for n=1:200w(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n); endfigure(1),plot(w);3、对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模型中间，分别记录两个时刻的波前快照（即区域内所有网格点的波场值）。
第一时刻为地震波还未传播到边界上的某时刻。
由稳定条件，设v为2000，可以令DT=0.001,DH=5,此时精度较高，且满足频散关系程序，图形如下：#include <stdio.h>#include <math.h>#include <stdlib.h>#define PI 3.141593#define FM 30#define R 3#define KN 200#define XN 101#define ZN 101#define DH 5#define DT 0.001void main(){FILE *fp;int i,j,k,m,n;float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],f[XN][ZN]; //不能直接初值为0 float u5[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;for(k=0;k<KN;k++)w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);for(i=0;i<XN;i++) //定义f函数，当且仅当i,j同时为50时，f为1，其余为0for(j=0;j<ZN;j++){if(i==50&&j==50)f[i][j]=1;else f[i][j]=0;}for(i=0;i<XN;i++)for(j=0;j<ZN;j++){ u1[i][j]=0.0;u2[i][j]=0.0;u3[i][j]=0.0;u4[i][j]=0.0;v[i][j]=2000; //速度相同表示同一介质}for(k=0;k<KN;k++){for(i=2;i<XN-2;i++)for(j=2;j<ZN-2;j++){uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];}for(m=0;m<XN;m++)for(n=0;n<ZN;n++){u1[m][n]=u2[m][n];u2[m][n]=u3[m][n];}if(k==100)for(m=0;m<XN;m++)for(n=0;n<ZN;n++)u4[m][n]=u3[m][n];//记录波前快照,中间点}if((fp=fopen("wavefront.dat","w"))!=NULL){fprintf(fp,"%d\n",XN);fprintf(fp,"%d\n",ZN);for(i=0;i<XN;i++)for(j=0;j<ZN;j++)fprintf(fp,"%f\n",u4[i][j]);fclose(fp);}}第二时刻为地震波已经传播到边界上的某时刻，体会其人工边界反射；程序图形如下：#include <stdio.h>#include <math.h>#include <stdlib.h>#define PI 3.141593#define FM 30#define R 3#define KN 200#define XN 101#define ZN 101#define DH 5#define DT 0.001void main(){FILE *fp;int i,j,k,m,n;float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],f[XN][ZN]; //不能直接初值为0 float u5[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;for(k=0;k<KN;k++)w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);for(i=0;i<XN;i++) //定义f函数，当且仅当i,j同时为50时，f为1，其余为0for(j=0;j<ZN;j++){if(i==50&&j==50)f[i][j]=1;else f[i][j]=0;}for(i=0;i<XN;i++)for(j=0;j<ZN;j++){ u1[i][j]=0.0;u2[i][j]=0.0;u3[i][j]=0.0;u4[i][j]=0.0;v[i][j]=2000; //速度相同表示同一介质}for(k=0;k<KN;k++){for(i=2;i<XN-2;i++)for(j=2;j<ZN-2;j++){uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];}for(m=0;m<XN;m++)for(n=0;n<ZN;n++){u1[m][n]=u2[m][n];u2[m][n]=u3[m][n];}if(k==160) //只需改动K值，即显示边界反射for(m=0;m<XN;m++)for(n=0;n<ZN;n++)u4[m][n]=u3[m][n];//记录波前快照，边界}if((fp=fopen("wavefront.dat","w"))!=NULL){fprintf(fp,"%d\n",XN);fprintf(fp,"%d\n",ZN);for(i=0;i<XN;i++)for(j=0;j<ZN;j++)fprintf(fp,"%f\n",u4[i][j]);fclose(fp);}}4、对于大模型，定义为水平层状速度模型；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的波前快照，体会地震波在两种介质的分界面上传播规律，指出哪是反射波，哪是透射波；这时取小模型的常量，为减少频散，速度v至少为2400程序图形如下：#include <stdio.h>#include <math.h>#include <stdlib.h>#define PI 3.141593#define FM 30#define R 3#define KN 200#define XN 200#define ZN 100#define DH 5#define DT 0.001void main(){FILE *fp;int i,j,k,m,n;float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],u5[XN][ZN]; //不能直接初值为0float f[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;for(k=0;k<KN;k++)w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);for(i=0;i<XN;i++)for(j=0;j<ZN;j++){ u1[i][j]=0.0;u2[i][j]=0.0;u3[i][j]=0.0;u4[i][j]=0.0;f[i][j]=0.0;if (j<=30)v[i][j]=2400; //第一层速度为2400elsev[i][j]=3000; //第二层速度为3000}f[100][10]=1; //定义f函数，在100，10为1for(k=0;k<KN;k++){for(i=2;i<XN-2;i++)for(j=2;j<ZN-2;j++){uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];u5[i][k]=u1[i][10]; //地震记录}for(m=0;m<XN;m++)for(n=0;n<ZN;n++){u1[m][n]=u2[m][n];u2[m][n]=u3[m][n];}if(k==100)for(m=0;m<XN;m++)for(n=0;n<ZN;n++)u4[m][n]=u3[m][n];//记录波前快照}if((fp=fopen("wavefront.dat","w"))!=NULL){fprintf(fp,"%d\n",XN);fprintf(fp,"%d\n",ZN);for(i=0;i<XN;i++)for(j=0;j<ZN;j++)fprintf(fp,"%f\n",u4[i][j]);fclose(fp);}}反射波透射波二是合成一个地震记录，即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分别对应哪些波？这时取小模型的常量，为减少频散，速度v至少为2400程序图像如下：#include <stdio.h>#include <math.h>#include <stdlib.h>#define PI 3.141593#define FM 30#define R 3#define KN 200#define XN 200#define ZN 200#define DH 5#define DT 0.001void main(){FILE *fp;int i,j,k,m,n;float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],u5[XN][ZN]; //不能直接初值为0 float f[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;for(k=0;k<KN;k++)w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);for(i=0;i<XN;i++)for(j=0;j<ZN;j++){ u1[i][j]=0.0;u2[i][j]=0.0;u3[i][j]=0.0;u4[i][j]=0.0;f[i][j]=0.0;if (j<=30)v[i][j]=2400; //第一层速度为2400elsev[i][j]=3100; //第二层速度为3100}f[100][10]=1; //定义f函数，在100，10为1for(k=0;k<KN;k++){for(i=2;i<XN-2;i++)for(j=2;j<ZN-2;j++){uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];u5[i][k]=u1[i][10]; //地震记录}for(m=0;m<XN;m++)for(n=0;n<ZN;n++){u1[m][n]=u2[m][n];u2[m][n]=u3[m][n];}if(k==100)for(m=0;m<XN;m++)for(n=0;n<ZN;n++)u4[m][n]=u3[m][n];//记录波前快照}if((fp=fopen("sei_record.dat","w"))!=NULL){fprintf(fp,"%d\n",XN/2);fprintf(fp,"%d\n",KN);for(i=0;i<XN;i=i+2)for(j=0;j<KN;j++)fprintf(fp,"%f\n",u5[i][j]);fclose(fp);} }从上向下看的两条直线为与震源同深度点的各点所有时刻的波场值，在向下与两直线紧邻的双曲线为第二层界面反射的记录。
