维波动方程第一类吸收边界条件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
(完整)吸收边界条件局域ABC
(3)保证算法稳定。
FDTD差分方程是不能用于截断边界的,因为它需要边界外面的点。根据实现ABC的途径,ABC分为如下三大类:
全域ABC
根据Green定理可知, 时刻边界上的场可以用 期间边界内部的所有点的场表示。这种ABC是严格的,但它涉及空间所有点在整个时间进程中的场值,即是“全域”的,所以,计算复杂,占内存大。实际模拟中很少采用。
(5—5b)
式中, .
一阶Taylor近似
(5-6)
适用于 非常小的情况. 非常小意味 非常小,即外向波沿 方向变化很小,等效地,波以基本上垂直于外边界的角度入射到外边界上.
(5—6)代入(5-5b),得
(5—7)
即一阶精度的单向波动方程
(5—8)
这正是沿 方向传播的波方程。
为了检验(5—8)的吸收特性,设
采用与前面类似的方法,可以得到反射系数为
同样,当 时反射为零, 越大反射越大, 时,全反射。
但二阶近似比一阶近似反射特性有所改善。
同理,我们可以得到三维情况的二阶单向波动方程
(5—12)
5.1。2 Mur有限差分方案
考虑(5—11),以 网格边界为例.波从 区域入射到 的边界,设 表示位于 处并与Yee网格边界相切的 或 直角坐标分量.Mur利用围绕网格点 展开的数值中心差分实现(5—11)中的偏导数的数值差分近似。
以直角坐标系中的二维波动方程为例
(5—2)
其中,偏微分算子L为
(5-3)
(5-2)的解为
式中 、 、 、 表示沿四个不同方向传输的波,以及
利用算子的因式分解,可以得到
(5—4)
式中, (5—5a)
波动方程一维波动方程的解
波动方程一维波动方程的解波动方程是描述物体在空间中传播波动的数学模型。
一维波动方程常用于描述沿直线传播的波动现象。
波动方程可以用以下形式表示:∂²u/∂t² = v² ∂²u/∂x²其中,u表示波动的位移或振幅,t表示时间,x表示空间位置,v 表示波速。
在解一维波动方程之前,我们先来讨论一下边界条件和初值条件,在实际问题中,这些条件通常会给出。
常见的条件有以下几种:1. 自由边界条件(自由端):在波动方程的一个或两个端点上,波动没有受到任何的约束或影响,即边界条件是自由的。
2. 固定边界条件(固定端):在波动方程的一个或两个端点上,波动被固定或限制住,不允许产生位移。
3. 开放边界条件:在波动方程的一个或两个端点上,波动可以自由流出或流入,即边界允许有反射和透射。
基于以上边界条件和初值条件,我们将根据不同场景进行求解一维波动方程的特定形式。
1. 矩形脉冲波动考虑一个在x轴上传播的矩形脉冲波动场景,即在初始时刻t=0处只有一个脉冲波动存在。
我们可以将初始条件表示为:u(x,0) = A (0<=x<=L), 其他地方 u(x,0) = 0其中,A表示波动的振幅,L表示脉冲波的长度。
解法:为了求解这个问题,我们可以使用方法之一-分离变量法。
首先,我们先猜测解的形式为u(x,t)=X(x)T(t),将其代入波动方程中,得到以下形式:X(x)T''(t) = v² X''(x)T(t)我们可以将该方程化简为两个简单的常微分方程,即:X''(x)/X(x) = T''(t)/(v²T(t)) = -λ²其中,λ是常数。
解得X(x) = Asin(λx) + Bcos(λx),T(t) = Csin(ωt) + Dcos(ωt)对于边界条件u(0,t) = u(L,t) = 0,可以得到X(0) = X(L) = 0。
波动方程
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)の解。
一维mur吸收边界条件fdtd形式
一维mur吸收边界条件fdtd形式一维Mur吸收边界条件在时域有限差分法(FDTD)中被广泛应用。
该边界条件用于模拟电磁波在计算域的边界处自由传播、消散或反射的情况。
本文将介绍一维Mur吸收边界条件的原理、应用以及其在FDTD中的形式。
一维Mur吸收边界条件是由J. E. Mur于1981年提出的,它基于完美匹配层(perfectly matched layer, PML)的思想。
PML是一种特殊的吸收层,能够有效地吸收电磁波并防止波的反射。
在FDTD中,Mur吸收边界条件通常用于近似计算域的边界情况,使得计算结果更加准确。
一维Mur吸收边界条件的形式如下:$$H_z(nx+1) = H_z(nx) - \frac{\Deltat}{\mu(n+\frac{1}{2})\Delta x}\left[ E_y(nx+1) - E_y(nx)\right]$$其中,$H_z(nx)$表示时刻n时空间节点nx处的磁场强度,$E_y(nx)$表示时刻n时空间节点nx处的电场强度,$\Delta t$为时间步长,$\mu$为磁导率,$\Delta x$为空间步长。
这种形式的一维Mur吸收边界条件通过将场量在计算域边界处的正常导数近似为差分形式来实现吸收。
它的核心思想是通过引入一系列的附加节点,对计算域的边界进行近似。
这些附加节点的值根据Mur 吸收边界条件的差分形式进行计算,从而消除了边界上的反射。
一维Mur吸收边界条件在FDTD中有着广泛的应用。
它可以用于模拟各种电磁波的传播过程,包括自由空间传播、波在介质中的传播以及波的反射和吸收等情况。
通过使用Mur吸收边界条件,可以有效地模拟出电磁波在计算域中的行为,并得到准确的计算结果。
然而,一维Mur吸收边界条件也有其局限性。
首先,它只适用于一维情况,不能直接推广到二维或三维情况。
其次,Mur吸收边界条件的效果受到许多参数的影响,如时间步长、空间步长、吸收层厚度等。
一维波动方程matlab
一维波动方程matlab一维波动方程是描述一维波动现象的数学模型。
在MATLAB中,我们可以使用偏微分方程求解器(Partial Differential Equation Toolbox)来求解一维波动方程。
假设我们有一个一维波动方程:∂^2u/∂t^2 = c^2 ∂^2u/∂x^2。
其中,u是波函数,t是时间,x是空间坐标,c是波速。
首先,我们需要定义波动方程的边界条件和初始条件。
然后,我们可以使用pdepe函数来求解一维波动方程。
以下是一个简单的MATLAB代码示例:matlab.function oned_wave_equation.x = linspace(0, 1, 100); % 定义空间网格。
t = linspace(0, 1, 100); % 定义时间网格。
m = 0; % 边界条件。
xmesh = x;tspan = t;sol =pdepe(m,@wave_eqn,@wave_initial,@wave_bc,xmesh,tspan); u = sol(:,:,1);surf(x,t,u);xlabel('Distance x');ylabel('Time t');zlabel('Wave function u');end.function [c,f,s] = wave_eqn(x,t,u,DuDx)。
c = 1; % 波速。
f = DuDx;s = 0;end.function u0 = wave_initial(x)。
u0 = sin(pix); % 初始波函数。
end.function [pl,ql,pr,qr] = wave_bc(xl,ul,xr,ur,t)。
pl = ul; % 左边界条件。
ql = 0;pr = 0; % 右边界条件。
qr = 1;end.在这个示例中,我们首先定义了空间和时间的网格,然后使用pdepe函数求解一维波动方程。
教材里面的附录 一维黎曼问题的求解编程代码 - 副本
计算得到的激波面很陡、很窄,计算激波精度是很高的。采用带开关函数的前置 人工滤波法能消除激波附近的非物理振荡,计算效果很好。 从图 A.2 和图 A.3 中可以看出通过激波后气体的密度、压力和速度都是增加 的;在压力分布中存在第二个台阶,表明在这里存在一个接触间断,在接触间断 两侧压力是有间断的, 而密度和速度是相等的。这个计算结果正确地反映了一维
图 A.1 激波管问题示意图
设气体是理想气体。一维 Riemann 问题在数学上可以用一维可压缩无黏气体
Euler 方程组来描述。在直角坐标系下量纲为一的一维 Euler 方程组为:
u f 0, 1 x 1 t x
其中 (A.1)
u u u , f u 2 p E ( E p)u
(A.4)
t 。计算实践表明, MacCormack 两步差分格式不能抑制激波附近非物 x
1 n n n n uin , j ui , j ui 1, j 2ui , j ui 1, j 2
理振荡。因此在计算激波时,必须采用人工黏性滤波方法: (A.5)
为了在激波附近人工黏性起作用,而在光滑区域人工黏性为零, 需要引入一个与 密度(或者压力)相关的开关函数:
MacCormack 两步差分格式:
uj u
其中 r
n 1 j n 1 2 n n un j r f j f j 1
1 1 n n 1 1 n 1 2 2 2 un u r f f j 1 j j j 2 2
(A.7)
其中参数 往往需要通过实际试算来确定,也可采 | 1 x x
(A.8)
一维波动方程的解法
一维波动方程的解法波动现象是自然界和人类生活中广泛存在的一种现象,它具有许多重要的物理意义,例如声波、光波等。
一维波动方程是描述波动的重要方程之一,本文将介绍一维波动方程的解法。
一、一维波动方程的基本形式和意义一维波动方程的基本形式为:$$\frac{\partial^2 u}{\partial t^2}-c^2\frac{\partial^2 u}{\partialx^2}=0$$其中,$u(x,t)$表示波动的幅度,$c$表示波速。
这个方程描述了介质中的一种波动现象:波动传播速度为$c$,波动在媒质中沿$x$轴方向的传播,波动的幅度随时间$t$的变化而变化。
在声波和电磁波中,$u$分别是空气压力和电场强度,$c$分别是声速和光速。
二、1. 分离变量法分离变量法是一种基本的解法,其思想是将波动方程中的未知函数$u(x,t)$表示成仅包含$x$的函数和仅包含$t$的函数的乘积形式:$$u(x,t)=X(x)T(t)$$将$u(x,t)$代入一维波动方程中,得到:$$\frac{T''(t)}{c^2T(t)}=\frac{X''(x)}{X(x)}=-\lambda^2$$其中,$\lambda$是一个常数。
由此可得到两个关于未知函数的简单微分方程:$$T''(t)+\lambda^2c^2T(t)=0$$和$$X''(x)+\lambda^2X(x)=0$$其中,第一个微分方程的解为:$$T(t)=A\cos(\lambda ct)+B\sin(\lambda ct)$$其中,A、B是常数。
第二个微分方程的解为:$$X(x)=C\cos(\lambda x)+D\sin(\lambda x)$$其中,C、D是常数。
因此,一维波动方程的通解为:$$u(x,t)=\sum_{n=1}^\infty a_n\cos(\lambda_n x)\cos(\lambda_n ct)+b_n\sin(\lambda_n x)\sin(\lambda_n ct)$$其中,$\lambda_n=n\pi/L$,$L$为介质的长度,$a_n$和$b_n$是待定常数。
MATLAB编辑一维波动方程的模拟的程序.doc
一维波动方程的模拟程序:function wave_equation() %一维线性齐次波动方程options={'空间杆长L','空间点数N' ,'时间点数M','波的相速度v',... '稳定条件的值r(取值必须小于1)','初始速度调用形式form(选择1或2)'}; topic='seting';lines=1;def={'1','100','100','1','1','1'};p=inputdlg(options,topic,lines,def);L=eval(p{1});N=eval(p{2});M=eval(p{3});v=eval(p{4});r=eval(p{5});%r的值必须小于1form=eval(p{6});%*************************************************************** h=L/N;%空间步长x=0:h:L;x=x';tao=r*h/v;%时间步长tm=M*tao;%波传播的总时间tmt=0:tao:tm;t=t';%计算边值和初值U=zeros(N+1,M+1);Uo=border_funo(t);Ue=border_fune(t);Ui=init_fun1(x);dUi=init_fun2(x);U(1,:)=Uo;U(N+1,:)=Ue;U(:,1)=Ui;if form==1U(:,2)=init_fun1(x)+tao*init_fun2(x);elsefor i=2:NU(i,2)=(1-r^2)*Ui(i)+0.5*r^2*(Ui(i+1)+Ui(i-1))+tao*dUi(i);endend%用差分法求解波动方程for j=3:(M+1)for i=2:NU(i,j)=2*(1-r^2)*U(i,j-1)+r^2*(U(i+1,j-1)+U(i-1,j-1))-U(i,j-2);endend%设置空间网格for i=1:N+1T(i,:)=t;endfor j=1:M+1X(:,j)=x;end%绘制出立体图形,即U-x-t图形figure(1)mesh(T,X,U)xlabel('t');ylabel('x');zlabel('U');title 一维波动方程的U-x-t图%绘制出平面图形,即U-x图形figure(2)for k=1:M+1plot(x,U(:,k))hold onendxlabel('x');ylabel('U');title 一维波动方程的U-x图function y=border_funo(t)%z=0处的边界条件y=1+t.*0;returnfunction y=border_fune(t)%z=L处的边界条件y=t.*0;returnfunction y=init_fun1(x)%初始位移条件w=3*pi;%w=pi;%w=2*pi;y=sin(w*x);returnfunction y=init_fun2(x)%初始速度条件y=x.*(1-x);return运行情况:按“run”运行时,弹出窗口将图框中的相关数据更改为:点击图框中的“OK”,在“command window”中输出结果为:。
维波动方程第一类吸收边界条件c++实现代码
#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_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][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_ol d[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*st ept-t0));//stepx1*source(m*stept+lowt)+u[m][Csize/2][0];stepx1*source(m*stept+lo wt)}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][Csize-2]-(1.0/sq rt(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 boundary u_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 boundary u_current[Csize-1][j]=u_old[Csize-1][j]+u_old[Csize-2][j]-u_past[Csize-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);}。
一种有效吸收边界条件的MATLAB实现
一种有效吸收边界条件的MATLAB 实现陈敬国中国地质大学(北京) 地球物理与信息技术学院 (100083)E-mail: chenjg_cugb@摘要:用有限差分法模拟地震波场是研究地震波在地球介质中传播的有效方法。
但我们在实验室进行波场数值模拟时有限差分网格是限制在人工边界里面,即引入了人工边界条件。
本文采用Clayton_Engquist_Majda 二阶吸收边界条件,通过MATLAB 编程实现了这一算法。
依靠MATLAB 具有更加直观的、符合大众思维习惯的代码,为用户提供了友好、简洁的程序开发环境,方便同行们交流。
关键词:有限差分法,地震波场,数值模拟,吸收边界条件,MATLAB1. 引言用有限差分法模拟地震波场是研究地震波在地球介质中传播的有效方法[1]。
但我们在实验室进行波场数值模拟时,只能在有限的空间进行,所以有限差分网格是限制在人工边界里面,即引入了人为的边界条件。
这种人为边界条件的引入将对有限区域内的波场值的计算带来严重影响,所以必须进行特殊的边界处理。
边界条件处理的好坏直接影响地震正演模拟的最终效果。
本文中我们采用Clayton_Engquist_Majda 二阶吸收边界条件[2]。
被称作是第四代计算机语言的MATLAB 语言,利用其丰富的函数资源把编程工作者从繁琐的程序代码中解放出来。
MATLAB 用更加直观的、符合大众思维习惯的代码,为用户提供了友好、简洁的程序开发环境。
本文介绍运用MATLAB 实现带有吸收边界条件的地震波场数值模拟方法和步骤,便于同行们交流,亦可用于本科地震理论的教学中,让学生们在程序演示中理解地震波的传播规律。
2. Clayton_Engquist_Majda 二阶吸收边界条件我们给定二维标量声波波动方程(含震源):2222221(,)P P f x z 2Px z v ∂∂∂++=∂∂∂t(1) 式中:是声波波场,是声波速度,P (,,)P x z t v (,)v x z (,)f x z 是震源。
- 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);}(注:素材和资料部分来自网络,供参考。