计算电磁学---有限差分法

第一章 有限差分法一元函数泰勒公式：设函数()f x 在0x 处的某邻域内具有1n +阶导数，则对该邻域异于0x 的任意点x ，在0x 与x 之间至少存在一点ξ，使得()20000000()()()()()()()()()2!!n n n f x f x f x f x f x x x x x x x R x n '''=+-+-+⋅⋅⋅+-+其中，(1)10()()()(1)!n n n f R x x x n ξ++=-+ 二元函数的泰勒公式：设函数(,)z f x y =在点00(,)x y 的某一邻域内连续且有直到1n +阶连续偏导数，00(,)x h y k ++为此邻域内任意点，则有00000020000100(,)(,)()(,)1()(,)2!1()(,)!1()(,)(1)!n n f x h y k f x y h k f x y x yh k f x y x y h k f x y n x y h k f x h y k n x yθθ+∂∂++=++∂∂∂∂+++⋅⋅⋅∂∂∂∂++∂∂∂∂+++++∂∂式中01θ<<；0000(,)0(,)mm m p p m p m x y p m p p f h k f x y c h kxy x y --=⎛⎫∂∂∂+= ⎪∂∂∂∂⎝⎭∑1.利用泰勒展开求不等间距的差分格式。

（1）2x y ϕ∂∂∂ （2）33xϕ∂∂解：(1)2600104001040101040010400010041()()2!11()()(,)!(1)!n n h h h h x y x yh h h h x h y h n x y n x yϕϕϕϕϕϕθθ+∂∂∂∂=+-++-++⋅⋅⋅∂∂∂∂∂∂∂∂+-++-+-+∂∂+∂∂（1.1）250020120020120102012002012000200121()()2!11()()(,)!(1)!n n h h h h x y x yh h h h x h y h n x y n x y ϕϕϕϕϕϕθθ+∂∂∂∂=+++++⋅⋅⋅∂∂∂∂∂∂∂∂++++++∂∂+∂∂（1.2）280090110090110109011009011000900111()()2!11()()(,)!(1)!n n h h h h x y x yh h h h x h y h n x y n x y ϕϕϕϕϕϕθθ+∂∂∂∂=+--+--+⋅⋅⋅∂∂∂∂∂∂∂∂+--+----∂∂+∂∂（1.3）270010030010030101003001003000100031()()2!11()()(,)!(1)!n n h h h h x y x yh h h h x h y h n x y n x yϕϕϕϕϕϕθθ+∂∂∂∂=+-+-+⋅⋅⋅∂∂∂∂∂∂∂∂+-+-+-∂∂+∂∂（1.4）()()()()605080702010400104020201200201202090110090110010*********()()2!1()()2!1()()2!1()(2!A B C D A h h A h h x y x y B h h B h h x y x y C h h C h h x y x y D h h D h h x y x ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ-+-+-+-∂∂∂∂=-++-+∂∂∂∂∂∂∂∂++++∂∂∂∂∂∂∂∂+--+--∂∂∂∂∂∂∂∂+-+-∂∂∂∂20)y ϕ （1.5）0102090100412011032222010209010222204012011030000Ah Bh Ch Dh Ah Bh Ch Dh Ah Bh Ch Dh Ah Bh Ch Dh -+-+=⎧⎪+--=⎪⎨+++=⎪⎪+++=⎩ （1.6） 有线性代数知识可知，有多种差分格式。

当090101002030110412h h h h h h h h =⎧⎪=⎪⎨=⎪⎪=⎩ 时，图1.1中左图退化为右图。

01020102040403032222010201022222040403030000Ah Bh Ch Dh Ah Bh Ch Dh Ah Bh Ch Dh Ah Bh Ch Dh -+-+=⎧⎪+--=⎪⎨+++=⎪⎪+++=⎩ （1.7）当点8在0点时，还可以退化为三个点。

图1.1(2) 33xϕ∂∂''3'2300100111()2!3!!nn n h h h R x n ϕϕϕϕϕϕ=++++⋅⋅⋅++''3'230200121212()()()()2!3!!nn n h h h h h h R x n ϕϕϕϕϕϕ=++++++⋅⋅⋅+++''3'2300300333()2!3!!nn n h h h R x n ϕϕϕϕϕϕ=++++⋅⋅⋅++由上述三式，可得：()()()()()102030'''233001013'''230012012'''233003033'222''3330112301123011()2611()2611()26()()222666A B C A h h B h h h h C h h A B C A B C A B C h h h h h h h h ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ-+-+-=+-⎡⎤+++-+⎢⎥⎣⎦+++⎡⎤⎡⎤=++++++-++-⎢⎥⎢⎥⎣⎦⎣⎦ 所以，2221123()0()0222A B C A B C h h h h ++=⎧⎪⎨+++=⎪⎩ 当C=1时，可求出A 、B 。

()()()()102030303331123666A B C A B C h h h h ϕϕϕϕϕϕϕ-+-+-=-++-图1.22. 等间距直接法（1）33x ϕ∂∂ （2）44xϕ∂∂解：（1） 由图1.3中的四点不能得出0点处中心差分格式，因此我们采用图1.4。

由0、1、2点得出1点的二阶导数差分格式，由0、3、4得出3点的二阶导数差分，再由1、3点的二阶导数差分格式得出0点的三阶导数差分格式。

具体推导如下：4330''''1''43301''''33104300122243123222222222h hh h h h h h h h hh ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ--⎛⎫- ⎪- ⎪⎛⎫-- ⎪ ⎪-⎝⎭⎝⎭===-+-+-=-+-=（2）由图1.4得：''''''''3001011243303001''''''''013033401300''''''300112433001''''2243300112343222233223h h h h h h h h h h h h h h h h ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ⎛⎫⎛⎫------ ⎪ ⎪- ⎪ ⎪⎛⎫⎛⎫-- ⎪ ⎪- ⎪ ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭⎝⎭===-+-+--+-+=--+== 30102134321432442h h h h h ϕϕϕϕϕϕϕϕϕϕ----+-++-=图 1.3图 1.43.二维正方形网格划分（中心差分情况）亥姆霍兹方程，得出其差分格式，并编写程序。

解：二维亥姆霍兹方程222220k x yϕϕϕ∂∂++=∂∂如图1.5所示，由21032202xx h ϕϕϕϕ⎛⎫-+∂= ⎪∂⎝⎭，22042202y y h ϕϕϕϕ⎛⎫-+∂= ⎪∂⎝⎭，得二维亥姆霍兹方程的0点处中心差分情况：2103204022220xyk h h ϕϕϕϕϕϕϕ-+-+++=当x y h h =时，上式变为：()1324022+4h k ϕϕϕϕϕ++=-图1.5程序实现：取k=0;采用SOR 迭代。

程序语言：VC++。

#include<iostream> #include<cmath> #include<cstdlib> #include<fstream> using namespace std;//迭代终止条件const double Ee0=0.00001;//误差限 const int NN=100000;//最大计算步数 const double w0=1.17;//超松弛因子//全局变量声明double a_l=0.04;//a 的长度 double b_l=a_l;//b 的长度 double dx=0.04/4;//空间间隔 double dy=dx;//常量声明int NX=int(a_l/dx)+1; int NY=int(b_l/dy)+1;//输出电势的坐标 int i_out=0.25*a_l/dx;int j_out=0.25*b_l/dy;class em{public:void initial();//初始化和边界条件int cacl(int nx_min,int nx_max,int ny_min,int ny_max,double **xy);//差分计算double **bxy;//电势};//构造函数动态分配所需的内存，并初始化变量，强制边界条件void em::initial(){//动态分配内存并通过循环初始化变量bxy=(double**)new char[sizeof(double*)*NX];int i,j;for(i=0; i<NX; i++){bxy[i]=new double[NY];//初始化各变量for(j=0; j<NY; j++)bxy[i][j]=0;}//强制边界条件for(i=0; i<NX; i++)for(j=0; j<NY; j++){if(i==0)bxy[i][j]=0;if(i==NX-1)bxy[i][j]=0;if(j==0)bxy[i][j]=0;if(j==NY-1)bxy[i][j]=100;}}//差分计算函数int em::cacl(int nx_min,int nx_max,int ny_min,int ny_max,double **xy){int sum=1;//迭代终止标志int sum0=1;double bxyn=0;//保存上一时刻的场值double e_xy;int i,j;for(i=nx_min;i<=nx_max;i++)for(j=ny_min;j<=ny_max;j++){bxyn=xy[i][j];xy[i][j]=0.25*(xy[i+1][j]+xy[i][j+1]+xy[i-1][j]+xy[i][j-1]);//差分方程xy[i][j]=bxyn+w0*(xy[i][j]-bxyn);e_xy=fabs(xy[i][j]-bxyn);if(e_xy<Ee0)sum0=1;else sum0=0;sum=sum&&sum0;}return sum;}int main(){ofstream outfile1;outfile1.open("u.dat",ios::out);int n;int check;em em1;em1.initial();for(n=0;n<NN;n++){check=em1.cacl(1,NX-2,1,NY-2,em1.bxy);if(check==1)break;}cout<<"误差限："<<Ee0<<endl;if(n==NN)cout<<"计算未终止"<<endl;cout<<"计算步数： "<<n<<endl;cout<<"坐标("<<i_out<<","<<j_out<<")电势为："<<em1.bxy[i_out][j_out]<<endl; outfile1<<em1.bxy[i_out][j_out]<<endl; outfile1.close(); return 0;}第二章 FDTD推导TE 情况下的FDTD 格式。