Диссертация (Задачи дифракции электромагнитных волн на системе произвольно расположенных тел и экранов), страница 13

P.778-783.65.Siegel K. M., Goodrich R. F., Weston V. H. Comments on far fienldscattering from bodies of revolution // Appl. Sci. Res. – 1959. Section B, vol 8. P.1-8.66.Smirnov Yu. G. Method of Volume Singular Integral Equation forDetermination of Permittivity of Dielectric Body in a Waveguide // Proceedingsof Progsess in Electromagnetics Research Symposium (PIERS 2009). – M., 2009.67.Smirnov Yu. G., Tsupak A. A.

Integro-differential Equations of theVector Problem of Electromagnetic Wave Diffraction by a System ofNonintersecting Screens and Inhomogeneous Bodies // Advances in MathematicalPhysics. -2015. –Vol. 2015. -Article ID 945965, 6 pages.68.Smirnov Yu. G., Tsupak A. A. Existence and uniqueness theorems inelectromagnetic diffraction on systems of lossless dielectrics and perfectlyconducting screens // Applicable Analysis. - 2016.

- Vol. 2016. - doi:10.1080/00036811.2016.1188289.69.Toyoda I., Matsuhara M., Kumagai N. Extended Integral EquationFormulation for Scattering Problems from a Cylindrical Scatterer // IEEETrans.Antennas Propagation. -1988. -Vol.36. -№ 11. P.1580-1586.70.Wang J. H. H. Generalized Moment Methods in Electromagnetics //New York: John Wiley & Sons, 1991.119Приложение 1В приложении представлен листинг основных блоков программы.Файл main.cpp#include <iostream>#include <fstream>#include "complex.h"#include "stream.h"#include "memory.h"#include "MatrixSystem.h"#include "Gradient.h"#include "FileOperations.h"using namespace std;// управляющие константыint num = 16; // кол-во конечных элементов сетки расположенных вдоль однойосиintN1 = 8*(num-1)*num + 4*num; // кол-во носителей в сетке на экранеintN2 = 3*(num-1)*(num)*(num); // кол-во носителей в сетке на телеintN = N1 + N2; //кол-во носителей в системе тело плюс экран плюспроволокаintip = 2; //количество точек интегририрования по стороне носителяdouble k = 1.0;double a = 0.0001;double Xi = 3.0*pi;double dim_s = 1.0;120double dim_b = 0.5;void GetGeometry(int num, int N, complex *W){for(int i=0; i<N; i++)W[i]=1.0;}void DivVector(int N, int N1, int N2, complex *U, complex *U1, complex *U2){int i,k=0,k1=0;for (i=0; i<N; i++){if (i<N1){U1[i] = U[i];}if ((i>=N1)&&(i<N1+N2)){U2[k] = U[i];k++;}}}void main(){//////////////////////cout<<"Start"<<endl;// выделяем память под матрицы и вектора..complex**A;CreateMatrixMemory(N, A);121complex *U;CreateVectorMemory(N, U, 0.0);complex *U1;CreateVectorMemory(N1, U1, 0.0);complex *U2;CreateVectorMemory(N2, U2, 0.0);complex *B;CreateVectorMemory(N, B, 0.0);complex *W;CreateVectorMemory(N, W, 1.0);GetGeometry(num, N, W);cout << "k = " << k << endl;cout << num << endl;CalculateMatrix(4, a, k, Xi, N1, N2, ip, num, dim_s, dim_b, A);//запись и чтение матрицы в бинарный файлcout<<"WriteMatrixToBinaryFile"<<endl;WriteMatrixToBinaryFile(num, N, A);cout<<"ReadMatrixToBinaryFile"<<endl;ReadMatrixToBinaryFile(num, N, A);cout<<"RightPart"<<endl;RightPart(k, N, N1, num, ip, dim_s, dim_b, B);// решаем СЛАУcout<<"Gradient"<<endl;122Gradient(N, A, B, W, U);DivVector(N, N1, N2, U, U1,U2);PrintVectorInExcelScreen(N1, num, U1);PrintVectorInExcelBody0Z( 0, num, U2);// удаляем памятьDeleteVectorMemory(N, W);DeleteVectorMemory(N, B);DeleteVectorMemory(N, U2);DeleteVectorMemory(N, U1);DeleteVectorMemory(N, U);DeleteMatrixMemory(N, A);cout<<"Finish"<<endl;///////////////////////}Файл Gradient.h#ifndef __GRADIENT_H#define __GRADIENT_H#include "complex.h"#include "memory.h"const double eps = 0.00000000000000000001;// квадрат погрешности// скалярное произведениеcomplex Norma2(int n, complex *x);123// умножение матрицы на векторvoid MatrixOnVector(int, int, complex **, complex *, complex *, complex *);// решение СЛАУ методом cопряженных градиентов для комплексной матрицыint Gradient(int, complex **, complex *, complex *, complex *);#endifФайл Gradient.cpp#include "gradient.h"#include "math.h"void MatrixOnVector(int n, int flag, complex **A, complex *B, complex *W,complex *U){int i,j;// умножение обычной матрицы на векторif (1 == flag){for(i=0; i<n; i++){U[i] = 0.0;for(j=0; j<n; j++){U[i] += A[i][j]*B[j]; // правая часть}}}else{// умножение сопряженной матрицы на векторfor(i=0; i<n; i++){U[i] = 0.0;for(j=0; j<n; j++){U[i] += conj(A[j][i])*B[j]; // правая часть124}}}for (i=0; i<n; i++){U[i] *= W[i];}}/* норма столбца */complex Norma2(int n, complex *x){int i;complex s,ret;complex *xs;xs = new complex[n];for (i=0; i<n; i++){xs[i] = real(x[i]) - _i*imag(x[i]);}s = 0.0;for (i=0; i<n; i++){s = s + x[i]*xs[i];}ret = s;125delete []xs;return ret;}/* метод сопряженных градиентов */int Gradient(int n, complex **A, complex *V, complex *W, complex *U){int i;//вспомагательные переменыеint k;// число итерацийcomplex b,s,al;//вспомагательные переменыеcomplex//вспомагательные переменыеh2,Nf02;//выделение памяти под вектораcomplex *f;CreateVectorMemory(n, f, 0.0);complex *f0;CreateVectorMemory(n, f0, 0.0);complex *w;CreateVectorMemory(n, w, 0.0);complex *q;CreateVectorMemory(n, q, 0.0);complex *h;CreateVectorMemory(n, h, 0.0);complex *p;CreateVectorMemory(n, p, 0.0);126// умножаем на вектор геометрииfor(i=0; i<n; i++){V[i] *= W[i];}//начальное приближениеfor(i=0; i<n; i++){f[i] = V[i];U[i] = V[i];}// правая часть нового матричного уравнения А1*Аu = А1*f и векторы w= A*u и q = A1*wMatrixOnVector(n, 0, A, f, W, f0);MatrixOnVector(n, 1, A, U, W, w);MatrixOnVector(n, 0, A, w, W, q);for(i=0; i<n; i++){h[i] = q[i] - f0[i];p[i] = h[i];}Nf02 = Norma2(n,f0);// скалярное произведениеh2 = Norma2(n, h);// скалярное произведениеk = 0;// запуск итерацийwhile (eps < abs(h2/Nf02)){127MatrixOnVector(n, 1, A, p, W, w);MatrixOnVector(n, 0, A, w, W, q);al = h2/Norma2(n,w);//новое приближение решенияfor(i=0; i<n; i++){U[i] = U[i] - al*p[i];h[i] = h[i] - al*q[i];}b = 1.0/h2;h2 = Norma2(n,h);b = b*h2;//новое значение pfor(i=0; i<n; i++){p[i] = h[i] + b*p[i];}cout << abs(h2/Nf02) << endl;k++;}cout << "Num iter = " << k << endl;// высвобождение памятиDeleteVectorMemory(n, p);DeleteVectorMemory(n, h);DeleteVectorMemory(n, q);128DeleteVectorMemory(n, w);DeleteVectorMemory(n, f0);DeleteVectorMemory(n, f);return 1;}Файл IntegralBodyScreen.h#include "complex.h"#include "special_fun.h"#include "grid_screen.h"#include "grid_body.h"complex Integral_BS2(double k, int ip, int num, double dim_s, double dim_b, int i,int j);complex Integral_BS3(double k, int ip, int num, double dim_s, double dim_b, int i,int j);complex Integral_Body_Screen(double k, int ip, int num, double dim_s, doubledim_b, int i, int j);Файл IntegralBodyScreen.cpp#include "IntegralBodyScreen.h"complex Integral_BS2(double k, int ip, int num, double dim_s, double dim_b, int i,int j){int i1, i2, i3;int j1, j2, j3;129int q1, q2, q3, q;q1 = GetNumber(1, j, num); q2 = GetNumber(2, j, num); q3 = GetNumber(3, j,num); q = GetNumber(4, j, num); // bodycomplex s, s1, s2, s3, s4;double h_s, len_s, h_b, len_b;//размеры сеткиlen_b = dim_b/((double)(num));h_b = len_b/((double)(ip));len_s = dim_s/((double)(num));h_s = len_s/((double)(ip));//screendouble *xp1; double *xp2; double *xp3; // массивы внутренних точек вячейки П+ для ребра рdouble *xm1; double *xm2; double *xm3; // массивы внутренних точек вячейки П- для ребра рxp1 = new double [ip];xm1 = new double [ip];xp2 = new double [ip];xm2 = new double [ip];xp3 = new double [ip];xm3 = new double [ip];// получаем массив внутренних точекGetInternalPointsS(i, num, ip, xp1,xp2,xp3, xm1,xm2,xm3); // screen//body130double *yp1; double *yp2; double *yp3; // массив внутренних точек вячейки П+ для ребра qdouble *ym1; double *ym2; double *ym3; // массив внутренних точек вячейки П- для ребра qyp1 = new double [ip]; ym1 = new double [ip];yp2 = new double [ip]; ym2 = new double [ip];yp3 = new double [ip]; ym3 = new double [ip];// получаем массив внутренних точекGetInternalPoints_Body(q1, q2, q3, q, dim_b, num, ip, yp1, yp2, yp3, ym1,ym2, ym3); // bodyint ap[3], am[3];TypeOfSupp(i, num, ap, am);ap[0] *= (-ip+1); ap[1] *= (-ip+1); ap[2] *= (-ip+1);am[0] *= (-ip+1); am[1] *= (-ip+1); am[2] *= (-ip+1);// интегрирование//cout<<endl<<endl<<"++";s1 = 0.0;for (j1=0; j1<ip; j1++){for (j2=0; j2<ip; j2++){for (j3=0; j3<ip; j3++){for (i1=0; i1<ip+ap[0]; i1++){for (i2=0; i2<ip+ap[1]; i2++){for (i3=0; i3<ip+ap[2]; i3++){131s1 += Baz_Body(num, dim_b, q1, q2, q3,q, yp1[j1], yp2[j2], yp3[j3])*Baz_p(0, i, num, xp1[i1], xp2[i2], xp3[i3], 4)*G(k,yp1[j1], yp2[j2], yp3[j3], xp1[i1], xp2[i2], xp3[i3]);}}}}}}//cout<<endl<<endl<<"-+";s2 = 0.0;for (j1=0; j1<ip; j1++){for (j2=0; j2<ip; j2++){for (j3=0; j3<ip; j3++){for (i1=0; i1<ip+ap[0]; i1++){for (i2=0; i2<ip+ap[1]; i2++){for (i3=0; i3<ip+ap[2]; i3++){s2 += Baz_Body(num, dim_b, q1, q2, q3,q, ym1[j1], ym2[j2], ym3[j3])*Baz_p(0, i, num, xp1[i1], xp2[i2], xp3[i3], 4)*G(k,ym1[j1], ym2[j2], ym3[j3], xp1[i1], xp2[i2], xp3[i3]);}}}}}132}//cout<<endl<<endl<<"+-";s3 = 0.0;for (j1=0; j1<ip; j1++){for (j2=0; j2<ip; j2++){for (j3=0; j3<ip; j3++){for (i1=0; i1<ip+am[0]; i1++){for (i2=0; i2<ip+am[1]; i2++){for (i3=0; i3<ip+am[2]; i3++){s3 += Baz_Body(num, dim_b, q1, q2, q3,q, yp1[j1], yp2[j2], yp3[j3])*Baz_p(1, i, num, xm1[i1], xm2[i2], xm3[i3], 4)*G(k,yp1[j1], yp2[j2], yp3[j3], xm1[i1], xm2[i2], xm3[i3]);}}}}}}//cout<<endl<<endl<<"--";s4 = 0.0;for (j1=0; j1<ip; j1++){for (j2=0; j2<ip; j2++){for (j3=0; j3<ip; j3++){133for (i1=0; i1<ip+am[0]; i1++){for (i2=0; i2<ip+am[1]; i2++){for (i3=0; i3<ip+am[2]; i3++){s4 += Baz_Body(num, dim_b, q1, q2, q3,q, ym1[j1], ym2[j2], ym3[j3])*Baz_p(1, i, num, xm1[i1], xm2[i2], xm3[i3], 4)*G(k,ym1[j1], ym2[j2], ym3[j3], xm1[i1], xm2[i2], xm3[i3]);}}}}}}s = s1+s2+s3+s4;delete []ym3;delete []yp3;delete []ym2;delete []yp2;delete []ym1;delete []yp1;delete []xm3;delete []xp3;134delete []xm2;delete []xp2;delete []xm1;delete []xp1;return s*h_b*h_b*h_b*h_s*h_s;}complex Integral_BS3(double k, int ip, int num, double dim_s, double dim_b, int i,int j){int i1, i2, i3;int j1, j2, j3;int q1, q2, q3, q;q1 = GetNumber(1, j, num); q2 = GetNumber(2, j, num); q3 = GetNumber(3, j,num); q = GetNumber(4, j, num); // Bodycomplex s, s1, s2, s3, s4;double h_s, len_s, h_b, len_b;//размеры сеткиlen_b = dim_b/((double)(num));h_b = len_b/((double)(ip));len_s = dim_s/((double)(num));h_s = len_s/((double)(ip));135//screendouble *xp1; double *xp2; double *xp3; // массивы внутренних точек вячейки П+ для ребра рdouble *xm1; double *xm2; double *xm3; // массивы внутренних точек вячейки П- для ребра рxp1 = new double [ip];xm1 = new double [ip];xp2 = new double [ip];xm2 = new double [ip];xp3 = new double [ip];xm3 = new double [ip];// получаем массив внутренних точекGetInternalPointsS(i, num, ip, xp1,xp2,xp3, xm1,xm2,xm3); // screen//bodydouble *yp1; double *yp2; double *yp3; // массив внутренних точек вячейки П+ для ребра qdouble *ym1; double *ym2; double *ym3; // массив внутренних точек вячейки П- для ребра qyp1 = new double [ip]; ym1 = new double [ip];yp2 = new double [ip]; ym2 = new double [ip];yp3 = new double [ip]; ym3 = new double [ip];// получаем массив внутренних точекGetInternalPoints_Body(q1, q2, q3, q, dim_b, num, ip, yp1, yp2, yp3, ym1,ym2, ym3);int ap[3], am[3];TypeOfSupp(i, num, ap, am);136ap[0] *= (-ip+1); ap[1] *= (-ip+1); ap[2] *= (-ip+1);am[0] *= (-ip+1); am[1] *= (-ip+1); am[2] *= (-ip+1);//cout<<endl<<endl<<"++";s1 = 0.0;for (j1=0; j1<ip; j1++){for (j2=0; j2<ip; j2++){for (j3=0; j3<ip; j3++){for (i1=0; i1<ip+ap[0]; i1++){for (i2=0; i2<ip+ap[1]; i2++){for (i3=0; i3<ip+ap[2]; i3++){s1 +=(-1.0/len_b)*(-1.0/len_s)*G(k,yp1[j1], yp2[j2], yp3[j3], xp1[i1], xp2[i2], xp3[i3]);}}}}}}//cout<<endl<<endl<<"+-";s2 = 0.0;for (j1=0; j1<ip; j1++){for (j2=0; j2<ip; j2++){for (j3=0; j3<ip; j3++){for (i1=0; i1<ip+ap[0]; i1++){for (i2=0; i2<ip+ap[1]; i2++){137for (i3=0; i3<ip+ap[2]; i3++){s2 += (1.0/len_b)*(-1.0/len_s)*G(k,ym1[j1], ym2[j2], ym3[j3], xp1[i1], xp2[i2], xp3[i3]);}}}}}}//cout<<endl<<endl<<"-+";s3 = 0.0;for (j1=0; j1<ip; j1++){for (j2=0; j2<ip; j2++){for (j3=0; j3<ip; j3++){for (i1=0; i1<ip+am[0]; i1++){for (i2=0; i2<ip+am[1]; i2++){for (i3=0; i3<ip+am[2]; i3++){s3 += (-1.0/len_b)*(1.0/len_s)*G(k,yp1[j1], yp2[j2], yp3[j3], xm1[i1], xm2[i2], xm3[i3]);}}}}138}}//cout<<endl<<endl<<"--";s4 = 0.0;for (j1=0; j1<ip; j1++){for (j2=0; j2<ip; j2++){for (j3=0; j3<ip; j3++){for (i1=0; i1<ip+am[0]; i1++){for (i2=0; i2<ip+am[1]; i2++){for (i3=0; i3<ip+am[2]; i3++){s4 += (1.0/len_b)*(1.0/len_s)*G(k,ym1[j1], ym2[j2], ym3[j3], xm1[i1], xm2[i2], xm3[i3]);}}}}}}s = s1+s2+s3+s4;delete []ym3;delete []yp3;delete []ym2;delete []yp2;139delete []ym1;delete []yp1;delete []xm3;delete []xp3;delete []xm2;delete []xp2;delete []xm1;delete []xp1;return s*h_b*h_b*h_b*h_s*h_s;}complex Integral_Body_Screen(double k, int ip, int num, double dim_s, doubledim_b, int i, int j){return -k*k*Integral_BS2(k, ip, num, dim_s, dim_b, i, j) + Integral_BS3(k, ip,num, dim_s, dim_b, i, j);}140Приложение 2В приложении представлена структура комплекса программ.Ввод: количество конечных элементов сетки, расположенных вдоль одной осиnum, тип integer; количество точек интегрирования по стороне носителя ip, тип integer; волновое число k, тип complex; размер тела dim_b, тип integer (рассматривается тело, имеющее формукуба со стороной dim_b); размер экрана dim_s, тип integer (рассматривается экран, имеющий формуцилиндра с основанием в виде квадрата, сторона квадрата - dim_s);Вывод: модули решения системы алгебраических уравнений.1.

Функции «CreateMatrixMemory» и «CreateVectorMemory» - выделениепамяти под матрицу и векторы соответственно.2. Функция «GetGeometry» - формирование вектора геометрии W,описывающего геометрию системы тел и экранов (если необходимо).3. Функция «CalculateMatrix» – формирование матрицы системы линейныхалгебраических уравнений; матрица состоит из четырех блоков, блокиматрицы заполняются подфункциями «Integral_Screen»,«Integral_Body_Screen», «Integral_Screen_Body», «Integral_Body»соответственно.4. Функция «RightPart» - заполнение вектора правой части СЛАУ.5.

Функция «Gradient» - решение СЛАУ итерационным методомсопряженных градиентов.6. Фукнция «PrintResults» - вывод полученных решений системы в файл«Results».141.