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

第一章 有限差分法 一元函数泰勒公式： 设函数在处的某邻域内具有阶导数，则对该邻域异于的任意点，在与之间至少存在一点，使得 其中， 二元函数的泰勒公式： 设函数在点的某一邻域内连续且有直到阶连续偏导数，为此邻域内任意点，则有 式中； 1.利用泰勒展开求不等间距的差分格式。 （1） （2） 解：(1) （1.1） （1.2） （1.3） （1.4） （1.5） （1.6） 有线性代数知识可知，有多种差分格式。 当 时，图1.1中左图退化为右图。 （1.7） 当点8在0点时，还可以退化为三个点。 图1.1 (2) 由上述三式，可得： 所以， 当C=1时，可求出A、B。 图1.2 2. 等间距直接法 （1） （2） 解： 由图1.3中的四点不能得出0点处中心差分格式，因此我们采用图1.4。 由0、1、2点得出1点的二阶导数差分格式，由0、3、4得出3点的二阶导数差分，再由1、3点的二阶导数差分格式得出0点的三阶导数差分格式。具体推导如下： （2）由图1.4得： 图 1.3 图 1.4 3.二维正方形网格划分（中心差分情况）亥姆霍兹方程，得出其差分格式，并编写程序。 解：二维亥姆霍兹方程 如图1.5所示， 由，，得二维亥姆霍兹方程的0点处中心差分情况： 当时，上式变为： 图1.5 程序实现：取k=0;采用SOR迭代。程序语言：VC++。 #includeiostream #includecmath #includecstdlib #includefstream 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; iNX; i++) { bxy[i]=new double[NY]; //初始化各变量 for(j=0; jNY; j++) bxy[i][j]=0; } //强制边界条件 for(i=0; iNX; i++) for(j=0; jNY; 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(x