维波动方程第一类吸收边界条件c++实现代码
教材里面的附录 一维黎曼问题的求解编程代码 - 副本
E[2]=(U[2]+p)*u; } /*------------------------------------------------------一维 MacCormack 差分格式求解器 入口: U， 上一时刻的 U 矢量，Uf、Ef，临时变量， dx，网格宽度，dt, 时间步长； 出口: U， 计算得到的当前时刻 U 矢量。 ---------------------------------------------------------*/ void MacCormack_1DSolver(double U[J+2][3],double Uf[J+2][3],double Ef[J+2][3],double dx,double dt) { int i,k; double r,nu,q; r=dt/dx; nu=0.25; for(i=1;i<=J;i++) { q=fabs(fabs(U[i+1][0]-U[i][0])-fabs(U[i][0]-U[i-1][0])) /(fabs(U[i+1][0]-U[i][0])+fabs(U[i][0]-U[i-1][0])+1e-100); //开关函数 for(k=0;k<3;k++) Ef[i][k]=U[i][k]+0.5*nu*q*(U[i+1][k]-2*U[i][k]+U[i-1][k]);//人工黏性项 } for(k=0;k<3;k++) for(i=1;i<=J;i++)U[i][k]=Ef[i][k]; for(i=0;i<=J+1;i++)U2E(U[i],Ef[i]); for(i=0;i<=J;i++) for(k=0;k<3;k++) Uf[i][k]=U[i][k]-r*(Ef[i+1][k]-Ef[i][k]); //U(n+1/2)(i+1/2) for(i=0;i<=J;i++)U2E(Uf[i],Ef[i]); //E(n+1/2)(i+1/2) for(i=1;i<=J;i++) for(k=0;k<3;k++) U[i][k]=0.5*(U[i][k]+Uf[i][k])-0.5*r*(Ef[i][k]-Ef[i-1][k]); //U(n+1)(i) } /*------------------------------------------------------输出结果, 用 Origin 数据格式画图
波动方程
1.1 波动方程的形式一维波动方程（描述弦的振动或波动现象的）()t x f x u a t u ,22222=∂∂-∂∂ 1.2 波动方程的定解条件（以一维波动方程为例）（1）边界条件①第一类边界条件（又称Dirichlet 边界条件）：弦振动问题中，弦的两端被固定在0=x 及l x =两点，因此有()0,0=t u ，()0,=t l u 。
②第二类边界条件（又称Neumann 边界条件）：弦的一端（例如0=x ）处于自由状态，即可以在垂直于x 轴的直线上自由滑动，未受到垂直方向的外力，此时成立0=∂∂=ox xu 。
也可以考虑更普遍的边界条件()t xu x μ=∂∂=0，其中()t μ是t 的已知函数。
③第三类边界条件：弦的一端固定在弹性支承上，不放考虑在l x =的一端，此时边界条件归结为0u =⎪⎭⎫ ⎝⎛+∂∂=l x u x σ。
也可以考虑更普遍的情况()t u x lx v u =⎪⎭⎫⎝⎛+∂∂=σ，其中()t v 是t 的已知函数。
1.3 利用叠加原理求解初值问题 初值问题()()()()⎪⎪⎩⎪⎪⎨⎧+∞<<∞=∂∂==+∞<<∞>=∂∂-∂∂)x -(,,:0t x 0,-t ,,22222x t u x u t x f x u a t u ψϕ （1） 利用叠加原理求解上述初值问题，叠加原理表明由()t x f ,所代表的外力因素和由()()x x ψϕ,所代表的初始振动状态对整个振动过程所产生的综合影响，可以分解为单独只考虑外力因素或只考虑初始振动状态对振动过程所产生的影响的叠加。
即如果函数()t x u ,1和()t x u ,2分别是下述初值问题（I ）()()()()⎪⎪⎩⎪⎪⎨⎧=∂∂===∂∂-∂∂2.1.....................,:0t 1.1. (022)222x t u x u x u a t u ψϕ和 （II ）()()()⎪⎪⎩⎪⎪⎨⎧=∂∂===∂∂-∂∂4.1....................................................................0,0:0t 3.1................................................................,22222t u u t x f x u a t u的解，那么()()t x u t x u u ,,21+=就一定是定解问题（1）的解。
一维波动方程
a a u 0 x t x t
ut au x v
按照第二章关于一阶线性偏微分方程的求解，我们也可以获得
D’Alembert公式(2.8).
§2 一维波动方程 7
《偏微分方程教程》 第四章 双曲型方程
上面对弦振动方程求解的特征线法, 亦适用于类似方程的
《偏微分方程教程》
第四章 双曲型方程
特征线法是求解一维双曲型方程Cauchy问题最基本的方法，
这个方法的实质是将方程沿特征线积分. 由第三章的特征概念知,
方程(2.1)的特征方程是
dx a dt 0
2 2 2
由此求得特征曲线为
c 其中 1 c2为任意常数. 为了将方程(2.1)化成第一标准型, 引入自变量变换
满足初始条件 u ( x 0) ( x)
ut ( x 0) ( x)
x x ,
(2.2)
其中a 是一个正常数，函数 ( x) C 2 ( x) C1 是定义在区 间
( ) 上的已知函数.
§2 一维波动方程 2
§2 一维波动方程
9
《偏微分方程教程》 第四章 双曲型方程
再对 求积分, 便得方程(2.11)的通解 u ( ) ( )d 1 ( ) 写成 其中1 ( )是 的任意函数. 若令 2 ( ) ( )d , 上式可
u( ) 1 ( ) 2 ()
其中 x 和 y 都是其变元的任意连续可微函数. 变回到原来的 变量1 和 2 , 便得到方程(2.10)的通解为
x u ( x y ) 1 ( xy ) xy 2 ( ) y
(2.14)
下面我们利用(2.10)中的初始条件来确定任意函数1 和 2.首先,容 易得到下面两个等式： 1 ( x) x2 ( x) ( x) (2.15)
一维传热问题边界条件的处理
一维传热问题边界条件处理当计算区域的边界为第二，第三类边界条件时，边界节点的温度是未知量。
为使内部节点的温度代数方程组得以封闭，有两类方法可以采用，即补充以边界节点代数方程的方法及附加原项法。
这里将介绍边界节点代数方程的方法。
对于无限大平板的第二类边界条件，采用泰勒展开法时，只要把边界条件B q x dX dT ==δλ中的导数用差分表达式来代替即可，即k q x T T B M M ⋅+=-δ111。
上式的截差为一阶，而内点上如采用中心差分，则截差为二阶。
为了得出具有二阶截差的公式，可以采用虚拟点法。
在边界外虚设一点M1+1，这样节点M1就可视为内节点，其一阶导数即可采用中心差分：B M M q xT T =--+δλ21111 为了消去TM1+1,由一维、稳态、含内热源的控制方程可得在M1点的离散形式：()02211111=++--+S x T T T M M M δλ从以上两式中消去11+M T 得，()()λδλδxq S x x T T B M M +∆+=-111其中2/x x δ=∆，是节点M1所代表的控制容积的厚度。
下面给出一个算例进行说明。
设有一导热型方程，022=-T dx T d ，边界条件为x=0,T=0; x=1, dT/dx=1。
试将该区域4等分，用区域离散方法求出各节点温度。
解：采用区域离散方法时，网格划分如下图所示，内点上采用中心差分。
右端点采用二阶截差，离散方程为： 0163332=-T T 01633432=-+-T T T 01633543=-+-T T T 41323354=+-T T编程解上述方程组得出每个节点的温度。
方程代码如下（Fortran6.6)：PROGRAM MAINUSE IMSLIMPLICIT NONEREAL :: A(4,4)=(/ 2.0625,-1.0,0.0,0.0,&-1.0,2.0625,-1.0,0.0,&0.0,-1.0,2.0625,-1.0,&0.0,0.0,-1.0,2.0625/) !矩阵A 的元素REAL :: B(4,1)=(/0.0,0.0,0.0,0.25/) !矩阵B 的元素REAL :: T(4,1) !4个节点的温度矩阵!EQUATION:!2.0625T2-T3=0!-T2+2.0625T3-T4=0!-T3+2.0625T4-T5=0!-T4+2.0625T5=0CALL LIN_SOL_GEN(A,B,T) !A*T=B,求解TWRITE(*,"(4F5.2)")TSTOPEND PROGRAM 0 T1 T3 T2 1/4 1/2 T5 T4 13/4。
一维波动方程的推导
一维波动方程的推导一维波动方程是描述一维介质中传播的波动现象的数学模型，它可以应用于声波、水波、电磁波等各种波动现象的研究。
其基本假设是介质中的波动是沿着介质传播的。
在推导一维波动方程时，我们需要先建立波动现象的数学模型。
假设介质中的波动是沿着x轴方向传播的，用u(x,t)表示波动处于x 点时的位移量。
我们需要考虑介质中的质点在时间t和t+Δt之间发生的位移量，即Δu(x,t)=u(x,t+Δt)-u(x,t)。
根据牛顿第二定律，质点在单位时间内所受到的合力等于质点的质量乘以加速度。
因此，介质中的质点在时间t和t+Δt之间的加速度可以表示为：a(x,t) = 1/ρ(x) * F(x,t)其中，ρ(x)是介质在x点处的密度，F(x,t)是介质在x点处的作用力。
根据胡克定律，介质中的质点在受到作用力时会发生弹性形变。
弹性形变的大小与作用力成正比，与介质的弹性系数成反比。
因此，介质在x点处的作用力可以表示为：F(x,t) = E(x) * u(x,t)/x其中，E(x)是介质在x点处的弹性系数，u(x,t)/x是介质在x点处的曲率。
将上述两个式子代入到a(x,t)的表达式中，得到：a(x,t) = 1/ρ(x) * E(x) * u(x,t)/x在介质中传播的波动是一种能量传输的过程。
波动在传播过程中，会带动介质中的质点振动，将能量从一个点传递到另一个点。
因此，介质中传播的波动在时间和空间上都是具有连续性的。
由此，我们可以得到波动方程的基本表达式：u(x,t)/t = c * u(x,t)/x其中，c=E/ρ，表示波动在介质中传播的速度的平方。
这就是一维波动方程的基本表达式。
在具体的应用中，我们需要根据不同的介质和波动特性，选择不同的初始条件和边界条件，来求解波动方程。
波动方程的吸收边界问题
波动方程的吸收边界问题波动方程是描述波动性现象重要的数学模型，涉及到横波与纵波的传播规律。
波动的传播不会受到边界的阻碍，因此，对于解决波动现象的数学模型中，吸收边界是一个非常重要的问题。
在实际应用中，吸收边界的概念是如下的：设计算区域外围为$U = {(x, y) | x\in [0, L],y\in [0, H]}$，则吸收边界是为了满足在$U$ 的封闭子集 $B = U \backslash I$ $(I\in [0, L]\times [0, H])$ 边界上的精确条件。
在波动方程的求解过程中，需要考虑对应于吸收边界的边界条件，以确保精确的计算结果。
现在介绍两种常用的吸收边界，分别是Mur吸收边界和Stefen 窄带边界条件。
1. Mur吸收边界Mur吸收边界是比较常见的一种吸收边界条件。
这种吸收边界的想法是模拟一种类似于黑洞的边界，能够吸收所有的波源、波浪和波波。
以二维波动方程为例，设波函数为 $u(x,y,t)$，则 Mur 条件中的 $x$方向边界为：$$u(x,y,t)=u(x_1-\Delta x,y,t)-R_x [u(x_1-\Delta x, y, t) - u(x, y, t)]$$其中，$\Delta x$ 为网格间距，$R_x$ 为吸收系数。
同理，$y$方向边界为：$$u(x,y,t)=u(x,y_1-\Delta y,t)-R_y [u(x, y_1-\Delta y, t) - u(x, y, t)]$$其中，$\Delta y$为网格间距，$R_x$ 为吸收系数。
Mur吸收边界的基本思想是在计算波函数时，将超过计算区域的波函数转化为一种相邻的波函数。
另外，使用 Mur 条件必须保证波函数的连续性，即在边界处存在连续性。
通过选定不同的吸收系数，可以控制边界对波函数的影响大小。
2. Stefen窄带边界条件Stefen窄带边界条件是另一种非常常用的吸收边界条件。
这种条件主要是通过对波函数进行变换，使得边界处的波函数能够逐渐减小，直至消失。
一维波动方程的解题方法及习题答案
第二篇数学物理方程—物理问题中的二阶线性偏微分方程及其解法 Abstracts:1、根据物理问题导出数理方程一偏微分方程； 2、给定数理方程的附加条件：初始条件、边界条件、物理条件 (自然条件，连接条件)，从而与数理方程一起构成定解问题； 3、方程齐次化； 4、数理方程的线性导致解的叠加。
一、数理方程的来源和分类(状态描述、变化规律)1、来源 I .质点力学：牛顿第二定律F =mr 连续体力学 II.麦克斯韦方程弹性体力学<(弹性定律)'弦 杆振动：出血力— a 2V 2 u (r , t ) = 0 (波动方程); 膜 0t 2 流体力学：质量守恒律:皿不V ・(p y ) = 0£ d t热力学物态方程：过+ (y -V )y ="p + f = 0 (Euler eq.).d t p JJ D .d c=fffp d i nV- D = p ; J E -d l =JJB -d s nVx E = B ;力B - d c= 0 nV- B = 0; J H - d l DjJ(j + D ) - d s nVx H = j + D . E = -V u , B = Vx A ,u ,A 满足波动方程。
、Lorenz 力公式n 力学方程;制axwell eqsT 电导定律n 电报方程。
IIL 热力学统计物理 热传导方程:以一 k V 2T = 0；特别：稳态(生= 0): V 23 = 0 (Laplace equation). < 八 01 扩散方程:0P - D V 2 p = 0. 、 01 IV.量子力学的薛定谔方程： i 方迦=—疟 V 2 u + Vu .0 01 2 m二、数理方程的导出推导泛定方程的原则性步骤：（1）定变量：找出表征物理过程的物理量作为未知数（特征量），并确定影响未知函数的自变量。
（2）立假设：抓主要因素，舍弃次要因素，将问题“理想化”--- “无理取闹”（物理趣乐）。
matlab练习程序（差分法解一维波动方程）
matlab练 习 程 序 （ 差 分 法 解 一 维 波 动 方 程 ）
实现了二维热传导方程数值解，这里我们计算波动方程数值解。 波动方程是一种双曲型偏微分方程。 这里依然用差分法计算。 一维波动方程如下：
写成差分形式：
整理一下就能得到u(i+1,j)。
en meshgrid(0:hx:x,0:ht:t); mesh(x1,t1,u)
结果如下：
matlab代码如下：
clear all;close all;clc;
t = 2;
%时间范围，计算到2秒
x = 1;
%空间范围，0-1米
m = 320;
%时间方向分320个格子
n = 64;
%空间方向分64个格子
ht = t/(m-1); %时间步长dt
hx = x/(n-1); %空间步长dx
u = zeros(m,n);
%设置边界条件 i=2:n-1; xx = (i-1)*x/(n-1); u(1,2:n-1) = sin(2*pi*xx); u(2,2:n-1) = sin(2*pi*xx);
%根据推导的差分公式计算 for i=2:m-1
for j=2:n-1 u(i+1,j) = ht^2*(u(i,j+1)+u(i,j-1)-2*u(i,j))/hx^2 + 2*u(i,j)-u(i-1,j);
波动方程的边值问题
波动方程的边值问题波动方程是物理学中非常重要的一个方程，它描述着物质振动的运动规律。
对于许多物理问题，我们经常需要解决它们的边值问题，而波动方程的边值问题也不例外。
首先，我们可以回忆一下波动方程的定义以及表达式。
波动方程可以描述波动物质（如弦、空气、水等）随时间的变化情况，它的一般形式为：$$\frac{\partial^2 u}{\partial t^2}=c^2\frac{\partial^2 u}{\partial x^2}$$其中 $u(x,t)$ 表示波动物质在位置 $x$ 和时间 $t$ 上的位移，$c$ 表示波速。
这个方程看起来很简单，但是实际上它涉及到许多复杂的物理学概念。
为了解决波动方程的边值问题，我们需要考虑它的初值条件和边界条件。
初值条件指的是在 $t=0$ 时刻波动物质的状态，而边界条件则是在波动物质的边界上的情况。
当我们给出了足够多的初值条件和边界条件，就可以通过一系列的计算来求解波动方程，并得到波动物质在不同时间和不同位置上的位移、速度等信息。
在波动方程的边值问题中，边界条件通常可以分为三类：第一类边界条件、第二类边界条件和第三类边界条件。
这些边界条件在实际问题中具有不同的意义和应用。
第一类边界条件（也称为 Dirichlet 边界条件）指定波动方程在边界上的取值，即 $u(x,0)=f(x)$ 和 $u(x,L)=g(x)$，其中 $L$ 是波动物质所在区域的长度。
这种边界条件意味着在边界上波动物质的位置是已知的，从而可以通过这些已知信息来解决波动方程的边值问题。
第二类边界条件（也称为 Neumann 边界条件）指定波动物质的导数在边界上的取值，即 $\frac{\partial u}{\partialx}(0,t)=h_1(t)$ 和 $\frac{\partial u}{\partial x}(L,t)=h_2(t)$。
这种边界条件通常用于描述波动物质的速度或者力的情况，它也可以用来解决一些特殊的问题。
《一维波动方程》课件
三维波动方程
描述空间波的传播
三维波动方程适用于描述在三维空间 中波的传播，例如声波、电磁波等。
物理应用广泛
三维波动方程在物理、工程等领域有 广泛的应用，如地震波传播、电磁波 传播等。
多维波动方程的解法
数值解法
对于多维波动方程，由于其复杂性， 通常采用数值解法来求解。常见的数 值解法包括有限差分法、有限元法等 。
解的物理意义
通过求解一维波动方程，可以得到波在空间中传播时的具体形式和性质，如波速、波长、振幅和相位 等。这些解具有明确的物理意义，可以用于描述和分析各种波动现象。
03
一维波动方程的解法
Chapter
分离变量法
总结词
通过将一维波动方程转化为多个常微分 方程，逐个求解，得到波动方程的解。
VS
详细描述
03
三维波动方程
描述三维空间中波的传播和变化规律，一般形式为：∂²u/∂t² = c² *
(∂²u/∂x² + ∂²u/∂y² + ∂²u/∂z²) + f(x, y, z, t)。
02
一维波动方程的建立
Chapter
一维波动方程的推导
波动现象的观察
波动现象在自然界中广泛存在，如声波、光波和水 波等。通过对这些现象的观察，可以发现波动具有 传播、干涉和衍射等特性。
有限差分法
总结词
通过将一维波动方程离散化，转化为差分方程组，然后求解差分方程组得到波 动方程的近似解。
详细描述
有限差分法是一种通过将一维波动方程离散化，转化为差分方程组的方法。在 离散化的过程中，需要考虑差分方程的稳定性和精度。然后利用数值计算方法 求解差分方程组，得到波动方程的近似解。
04
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
#include "stdafx.h"#include <iostream>#include <fstream>#include <cmath>#include<string>using namespace std;const double pi=4*atan(1.0);double freq=45;double sb=7.45;double t1=2*pi/(sb*4);double source(double t){//double t2=0.0;if(t<=t1) return (sin(sb*4*t-pi/2)+1)/10;else{ double tep=0.0; return tep;}//return((1-2*pi*pi*freq*freq*t*t)*exp(-pi*pi*freq*freq*t*t)+1);//Ricker子波}void update_Vn(double upt,double lowt,double upx1,double lowx1){int i,j,m;const int Csize=300;double deg=0;double stepx1=abs(upx1-lowx1)/(Csize-1);//double te=sqrt(static_cast<double>(3.0/8.0));double stept=sqrt(static_cast<double>(1.0/2.0))*stepx1/2.0;//int tn=static_cast<int>(upt/stept);double r=stept/stepx1;double **u_current,**u_old,**u_past;u_current=new double *[Csize];u_old=new double*[Csize];u_past=new double*[Csize];for(i=0;i<Csize;i++){u_current[i]=new double [Csize];u_old[i]=new double[Csize];u_past[i]=new double[Csize];}for(i=0;i<Csize;i++)for(j=0;j<Csize;j++){u_current[i][j]=0;u_old[i][j]=0;u_past[i][j]=0;}double ck[Csize][Csize];int flag=0;for(j=0;j<Csize;j++){for(i=0;i<Csize;i++){if(j<i) ck[i][j]=4;else ck[i][j]=1;}}}string str;cout<<"\n 输入保存文件名:";cin>>str;ofstream fout(str);if(!fout){cout<<"\n 不能打开文件"<<str<<endl;exit(1);}m=0;double f0=2.0/(stept*30);double t0=4.0/f0;while(m<1500&&((m++)*stept+lowt)<upt){for(i=0;i<Csize;i++)for(j=0;j<Csize;j++){if((i!=0&&(i!=Csize-1))&&(j!=0&&j!=(Csize-1)))//cope with the internal points of the interesting domain{//if(i==(Csize/2)&&j==(Csize/2))u_current[i][j]=(r*r/ck[i][j])*(u_old[i+1][j]+u_old[i][j+1]+u_old[i-1] [j]+u_old[i][j-1])-2*u_old[i][j]*(2*r*r/ck[i][j]-1.0)-u_past[i][j]//+stept*stept*source(m*stept+lowt);//stept*stept*exp(-f0*f0*(m*stept-t0)*(m*stept-t0));//stept*stept*source(m*stept+lowt);//elseu_current[i][j]=(r*r/ck[i][j])*(u_old[i+1][j]+u_old[i][j+1]+u_ol d[i-1][j]+u_old[i][j-1])-2*u_old[i][j]*(2*r*r/ck[i][j]-1.0)-u_past[i][ j];//u[m][i][j]=r*r*ck[i][j]*ck[i][j]*(u[m-1][i+1][j]+u[m-1][i][j+1]+u [m-1][i][j-1]+u[m-1][i-1][j])-2*(2*r*r*ck[i][j]*ck[i][j]-1)*u[m-1][i][ j]//-u[m-2][i][j];//u_current[i][j]=r*r*ck[i][j]*ck[i][j]*(u_old[i+1][j]+u_old[i][j+1]+u _old[i][j-1]+u_old[i-1][j])-2*(2*r*r*ck[i][j]*ck[i][j]-1)*u_old[i][j]//-u_past[i][j];//u_current[i][j]=(power(r,2.0)/ck[i][j])*(u_old[i+1][j]+u_old[i][j +1]+u_old[i-1][j]+u_old[i][j-1])-2*u_old[i][j]*(2*pow(r,2.0)/ck[[i][j] -1)-u_past[i][j];u_current[i][j]=(r*r/ck[i][j])*(u_old[i+1][j]+u_old[i][j+1]+u_old[i-1] [j]+u_old[i][j-1])-2*u_old[i][j]*(2*r*r/ck[i][j]-1.0)-u_past[i][j];}//u[m][0][j]=u[m][1][j];//u[m][Csize-1][j]=u[m][Csize-2][j];////u[m][i][0]=u[m][i][1];// u[m][i][Csize-1]=u[m][i][Csize-2];if(i==Csize/2){u_current[Csize/2][1]=u_current[Csize/2][0]+stepx1*exp(-f0*f0*(m*stept -t0)*(m*stept-t0));//stepx1*source(m*stept+lowt)+u[m][Csize/2][0];step x1*source(m*stept+lowt)}elseu_current[i][0]=u_old[i][0]+u_old[i][1]-u_past[i][1]+(1.0/sqrt(ck[i][0 ]))*r*(u_old[i][1]-u_old[i][0]-u_past[i][2]+u_past[i][1]);//top side absorobing boundary u[m][i][0]=u[m][i][1];//u_current[i][0]=u_old[i][0]+u_old[i][1]-u_past[i][1]+ck[i][0]*r*( u_old[i][1]-u_old[i][0]-u_past[i][2]+u_past[i][1]);//top side absorobing boundary u[m][i][0]=u[m][i][1];//u_current[i][Csize-1]=u_old[i][Csize-1]+u_old[i][Csize-2]-u_past[i][Cs ize-2]-(1.0/sqrt(ck[i][Csize-1]))*r*(u_old[i][Csize-1]-u_old[i][Csize-2]-u_past[i][Csize-2]+ u_past[i][Csize-3]);// down side absorbing boundaryu_current[0][j]=u_old[0][j]+u_old[1][j]-u_past[1][j]+(1.0/sqrt(ck[i ][Csize-1]))*r*(u_old[1][j]-u_old[0][j]-u_past[2][j]+u_past[1][j]);//Left side absorbing boundaryu_current[Csize-1][j]=u_old[Csize-1][j]+u_old[Csize-2][j]-u_past[Cs ize-2][j]-(1.0/sqrt(ck[i][Csize-1]))*r*(u_old[Csize-1][j]-u_old[Csize-2][j]-u_past[Csize-2][j]+u_past[Csize-3][j]);//Right side absorbing boundary}if(m%10==0){for(i=0;i<Csize;i++){for(j=0;j<Csize;j++){fout<<" "<<u_current[j][i];//fout<<" "<<ck[i][j];}fout<<endl;}fout<<endl;}for(i=0;i<Csize;i++)for(j=0;j<Csize;j++){u_past[i][j]=u_old[i][j];u_old[i][j]=u_current[i][j];u_current[i][j]=0;}// m++;}fout.close();for(i=0;i<Csize;i++) {delete[] u_current[i],u_old[i],u_past[i];}delete [] u_current,u_old,u_past;}void main () //主函数{double upt=17.8*t1*10;double lowt=0,upx1=3.0,lowx1=-3.0;// const int Csize=100;// const int tn=100;//double upt=1;update_Vn(upt, lowt,upx1,lowx1);}（注：素材和资料部分来自网络，供参考。