维波动方程第一类吸收边界条件c++实现代码
- 格式:doc
- 大小:25.50 KB
- 文档页数:7
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)的解。
一维传热问题边界条件处理当计算区域的边界为第二,第三类边界条件时,边界节点的温度是未知量。
为使内部节点的温度代数方程组得以封闭,有两类方法可以采用,即补充以边界节点代数方程的方法及附加原项法。
这里将介绍边界节点代数方程的方法。
对于无限大平板的第二类边界条件,采用泰勒展开法时,只要把边界条件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)立假设:抓主要因素,舍弃次要因素,将问题“理想化”--- “无理取闹”(物理趣乐)。
波动方程的边值问题波动方程是物理学中非常重要的一个方程,它描述着物质振动的运动规律。
对于许多物理问题,我们经常需要解决它们的边值问题,而波动方程的边值问题也不例外。
首先,我们可以回忆一下波动方程的定义以及表达式。
波动方程可以描述波动物质(如弦、空气、水等)随时间的变化情况,它的一般形式为:$$\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)$。
这种边界条件通常用于描述波动物质的速度或者力的情况,它也可以用来解决一些特殊的问题。
#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);}(注:素材和资料部分来自网络,供参考。