Програма для рішення інтегрального рівняння Фредгольма другого роду простим методом, заснованим на одночленних рекурентних формулах

Програма написана мовою програмування "Фортран" відповідно до алгоритму, приведеним у дипломній роботі.

Відзначимо деякі особливості, що були використані при написанні програм, приведених у додатках.

Уведення вихідних даних виробляється з файлу input.txt. Висновок результатів виробляється у файл result.txt.

Для наближеного обчислення визначених інтегралів виду

використовувався метод Симпсона (метод парабол) чисельного інтегрування [22,23] відповідно до якого:

,

де , 2m - число крапок розподілу інтервалу інтегрування.

Функція Ханкеля другого роду нульового порядку, рівна (x) = J ₀ {x)-i ₀ (x) представлялася через апроксимаційні багаточлени для функцій Бесселя і Неймана, що мають наступний вид [23] при 0 < х ≤ 3:

J₀(x) = l-2,2499997(x/3)² + l,2656208(x/3)⁴ -0,3163866(x/3)⁶ +

+ 0,0444479(x/3)⁸ -0,0039444(x/3)¹⁰ + 0,00021(x/3)¹²,

(x)=(2/π)ln(x/2) J _o(x)+0,36746691+0,60559366(x/3)²-,74350384(x/3)⁴+

+0,25300117(x/3)⁶-0,04261214(x/3)⁸+0,00427916(x/3)¹⁰ -0,00024846(x/3)¹²

 

при 3 ≤ х < ∞:

 

 

=0,79788456-0,00000077(3/ x)- 0,0055274(3/ x)² -0,00009512(3/ х)³-

-0,00054125(3/ x)⁴-0,00029333(3/x)⁵+0,00013558(3/x)⁶,

θ₀=х-0,78539816-0,04166397(3/х)-0,00003954(3/х)²+0,00262573(3/х) ³ -

-0,00054125(3 /x)⁴-0,00029333(3/x)⁵+0,00013558(3/x)⁶.

* Програма обчислення коефіцієнта придушення однорідної і неоднорідної структури

 

* визначення типу перемінних

character pop,pap,pup,pep,pok,pro,prr,ptt,prp,ppr,ptr

double precision yn,yk,pi,ba,t,hpd,shr,hpp,sg,ni,ki,an

* підготовка файлів для введення (уведення данних.txt) і висновку даних (результат.txt)

open(2,file='данные.txt',status='old')

open(1,file='etta_rez.xlt',status='new')

open(3,file='lambda.xlt',status='new')

open(4,file='impedance.xlt',status='new')

open(5,file='modul_J.xlt',status='new')

open(6,file='real_J.xlt',status='new')

open(7,file='im_J.xlt',status='new')

read(2,*)pop

read(2,*)nn

read(2,*)pap

read(2,*)yn

read(2,*)pup

read(2,*)yk

read(2,*)pep

read(2,*)ba

read(2,*)prr

read(2,*)hpp

read(2,*)pok

read(2,*)lp

read(2,*)pro

read(2,*)ni

read(2,*)ptr

read(2,*)ki

read(2,*)prp

read(2,*)sg

read(2,*)ptt

read(2,*)shr

read(2,*)ppr

read(2,*)hpd

pi=3.14159265359D0

t=(yk-yn)/nn

* Передача даних підпрограмі рішення інтегрального рівняння

 

call integr(nn,yn,yk,t,pi,ba,lp,hpd,shr,hpp,sg,ni,ki)

stop

end

 

* Підпрограма рішення інтегрального рівняння

* Визначення типу перемінних

subroutine integr(nn,yn,yk,t,pi,ba,lp,hpd,shr,hpp,sg,ni,ki)

complex*16 mat,f,fx,ro,jot1,dr,fy,u1,fun,fn,rez,han,han1,imp,ss,u2

double precision a0,a1,a2,a3,a4,a5,a6,b0,b1,b2,b3,b4,b5,b6,pi,kof

double precision d0,d1,d2,d3,d4,d5,z0,z1,z2,z3,z4,z5,z6,ba,t,ck

double precision am,h,a,x,yhn,yhk,at,b,d,q,p,ot,al,y1,yn,yk,ww,hpd

double precision shr,hpp,aa,ab,ac,ad,af,yy,xx,yx,fu,ff,y,yyy

double precision vvs,vas,vvas,sks,sg,ni,ki,bda

dimension mat(1000),jot1(1000,1000),ro(1000),fy(1000),ss(1000)

dimension u1(1000),fun(1000,1000),at(1000),ot(1000),fu(1000)

dimension han(1000,1000),han1(1000,1000),ww(1000),u2(1000)

dimension ff(1000),fn(1000),yyy(1000),rez(1000)

 

 

if(hpp-ba)22,23,24

22 write(1,*)'Неправильні дані кінцевого значення довжини хвилі'

return

24 kkk=(hpp-ba)/0.01d0

do 333 llll=1,kkk+1

if(hpp-ba)26,27,26

27 return

26 bda=ba+(llll-1)*0.01d0

goto 25

23 bda=ba

 

25 ck=2.d0*pi/bda

lll=((ki-ni)+sg)/sg

do 123 ll=1,lll

m=500

am=m

h=t/(2.d0*am)

m1=m-1

 

* Визначення коефіцієнтів для функцій Бесселя

a0=0.79788456d0

a1=-0.00000077d0

a2=-0.00552740d0

a3=-0.00009512d0

a4=0.00137237d0

a5=-0.00072805d0

a6=0.00014476d0

b0=-0.78539816d0

b1=-0.04166397d0

b2=-0.00003954d0

b3=0.00262573d0

b4=-0.00054125d0

b5=-0.00029333d0

b6=0.00013558d0

d0=-2.2499997d0

d1=1.2656208d0

d2=-0.3163866d0

d3=0.0444479d0

d4=-0.0039444d0

d5=0.00021d0

z0=0.36746691d0

z1=0.60559366d0

z2=-0.74350384d0

z3=0.25300117d0

z4=-0.04261214d0

z5=0.00427916d0

z6=-0.00024846d0

f=dcmplx(0.d0,1.d0)

sks=1.d0/(120.d0*bda)

* Блок визначення проходження елементарних проміжків

y1=yn+t/2.d0

do 1 j=1,nn

do 1300 i=1,nn

1300 yyy(i)=yn+t*(i-1)

yhn=yyy(j)

yhk=yhn+t

* Інтегрування функції Ганкеля по формулі Симпсона

* GANKEL-A

a=ck*(y1-yhn)

aa=ck*(y1-yhk)

if(a)50,51,51

50 a=dabs(ck*(yhn-y1))

aa=dabs(ck*(y1-yhk))

ab=ck*(yhk-y1)

ac=ab

ad=ck*(y1-yhn)

af=ck*(yhn-y1)

goto 52

51 if(aa)70,70,71

70 a=dabs(ck*(y1-yhn))

aa=dabs(ck*(y1-yhk))

ab=ck*(yhk-y1)

ac=ab

ad=ck*(y1-yhn)

af=ad

goto 52

71 a=dabs(ck*(y1-yhn))

aa=dabs(ck*(y1-yhk))

ab=ck*(yhk-y1)

ac=ck*(y1-yhk)

ad=ck*(y1-yhn)

af=ad

52 if(a -3.d0)4,4,5

4 if(aa-3.d0)72,72,5

72 x=af

yx=ac

xx=ab

yy=ad

at(j)=(d0/(3.d0**2))*((xx**3+yy**3)/3.d0)

at(j)=at(j)+(d1/(3.d0**4))*((xx**5+yy**5)/5.d0)

at(j)=at(j)+(d2/(3.d0**6))*((xx**7+yy**7)/7.d0)

at(j)=at(j)+(d3/(3.d0**8))*((xx**9+yy**9)/9.d0)

at(j)=at(j)+(d4/(3.d0**10))*((xx**11+yy**11)/11.d0)

at(j)=at(j)+(d5/(3.d0**12))*((xx**13+yy**13)/13.d0)

at(j)=(yhk-yhn)+at(j)/ck

ot(j)=(xx*(dlog(yx/2.d0)-1.)+yy*(dlog(x/2.d0)-1.d0))/ck

kof=(2.*d0*((2.d0/3.d0)**2))/(3.d0*ck)

ot(j)=ot(j)+(dlog(yx/2.d0)-1.d0/3.d0)*((xx/2.d0)**3)*kof

ot(j)=ot(j)+(dlog(x/2.d0)-1.d0/3.d0)*((yy/2.d0)**3)*kof

kof=(2.d0*d1*((2.d0/3.d0)**4))/(5.d0*ck)

ot(j)=ot(j)+(dlog(yx/2.d0)-1.d0/5.d0)*((xx/2.d0)**5)*kof

ot(j)=ot(j)+(dlog(x/2.d0)-1.d0/5.d0)*((yy/2.d0)**5)*kof

kof=(2.d0*d2*((2.d0/3.d0)**6))/(7.d0*ck)

ot(j)=ot(j)+(dlog(yx/2.d0)-1.d0/7.d0)*((xx/2.d0)**7)*kof

ot(j)=ot(j)+(dlog(x/2.d0)-1.d0/7.d0)*((yy/2.d0)**7)*kof

kof=(2.d0*d3*((2.d0/3.d0)**8))/(9.d0*ck)

ot(j)=ot(j)+(dlog(yx/2.d0)-1.d0/9.d0)*((xx/2.d0)**9)*kof

ot(j)=ot(j)+(dlog(x/2.d0)-1.d0/9.d0)*((yy/2.d0)**9)*kof

kof=(2.d0*d4*((2.d0/3.d0)**10))/(11.d0*ck)

ot(j)=ot(j)+(dlog(yx/2.d0)-1.d0/11.d0)*((xx/2.d0)**11)*kof

ot(j)=ot(j)+(dlog(x/2.d0)-1.d0/11.d0)*((yy/2.d0)**11)*kof

kof=(2.d0*d5*((2.D0/3.d0)**12))/(13.d0*ck)

ot(j)=ot(j)+(dlog(yx/2.d0)-1.d0/13.d0)*((xx/2.d0)**13)*kof

ot(j)=ot(j)+(dlog(x/2.d0)-1.d0/13.d0)*((yy/2.d0)**13)*kof

ot(j)=(2.d0/pi)*ot(j)

ot(j)=ot(j)+z0*(yhk-yhn)

ot(j)=ot(j)+(z1/((3.d0**2)*ck))*((xx**3+yy**3)/3.d0)

ot(j)=ot(j)+(z2/((3.d0**4)*ck))*((xx**5+yy**5)/5.d0)

ot(j)=ot(j)+(z3/((3.d0**6)*ck))*((xx**7+yy**7)/7.d0)

ot(j)=ot(j)+(z4/((3.d0**8)*ck))*((xx**9+yy**9)/9.d0)

ot(j)=ot(j)+(z5/((3.d0**10)*ck))*((xx**11+yy**11)/11.d0)

ot(j)=ot(j)+(z6/((3.d0**12)*ck))*((xx**13+yy**13)/13.d0)

fn(j)=dcmplx(at(j),-ot(j))

goto 1

5 a=dabs(ck*(yhn-y1))

x=3.d0/a

at(j)=a0+a1*x+a2*(x**2)+a3*(x**3)+a4*(x**4)+a5*(x**5)

at(j)=at(j)+a6*(x**6)

ww(j)=at(j)/dsqrt(a)

ot(j)=a+b0+b1*x+b2*(x**2)+b3*(x**3)+b4*(x**4)+b5*(x**5)

ot(j)=ot(j)+b6*(x**6)

at(j)=ww(j)*dcos(ot(j))

ot(j)=ww(j)*dsin(ot(j))

fu(j)=at(j)

ff(j)=ot(j)

* GANKEL-B

b=dabs(ck*(yhk-y1))

x=3.d0/b

at(j)=a0+a1*x+a2*(x**2)+a3*(x**3)+a4*(x**4)+a5*(x**5)

at(j)=at(j)+a6*(x**6)

ww(j)=at(j)/dsqrt(b)

ot(j)=b+b0+b1*x+b2*(x**2)+b3*(x**3)+b4*(x**4)+b5*(x**5)

ot(j)=ot(j)+b6*(x**6)

at(j)=ww(j)*dcos(ot(j))

ot(j)=ww(j)*dsin(ot(j))

fu(j)=fu(j)+at(j)

ff(j)=ff(j)+ot(j)

* GANKEL-(B-H)

d=dabs(ck*((yhk-h)-y1))

x=3.d0/d

at(j)=a0+a1*x+a2*(x**2)+a3*(x**3)+a4*(x**4)+a5*(x**5)+a6*(x**6)

ww(j)=at(j)/dsqrt(d)

ot(j)=d+b0+b1*x+b2*(x**2)+b3*(x**3)+b4*(x**4)+b5*(x**5)

ot(j)=ot(j)+b6*(x**6)

at(j)=ww(j)*dcos(ot(j))

ot(j)=ww(j)*dsin(ot(j))

fu(j)=fu(j)+4.d0*at(j)

ff(j)=ff(j)+4.d0*ot(j)

 

do 40 l=1,m1

al=l

* GANKEL-(A+(2*AK-1)*H)

q=dabs(ck*(yhn+(2.d0*al-1.d0)*h-y1))

x=3.d0/q

at(j)=a0+a1*x+a2*(x**2)+a3*(x**3)+a4*(x**4)+a5*(x**5)+a6*(x**6)

ww(j)=at(j)/dsqrt(q)

ot(j)=q+b0+b1*x+b2*(x**2)+b3*(x**3)+b4*(x**4)+b5*(x**5)

ot(J)=ot(j)+b6*(x**6)

at(j)=ww(j)*dcos(ot(j))

ot(j)=ww(j)*dsin(ot(j))

fu(j)=fu(j)+4.d0*at(j)

ff(j)=ff(j)+4.d0*ot(j)

* GANKEL-(A+2*AK*H)

p=dabs(ck*(yhn+2.d0*al*h-y1))

x=3.d0/p

at(j)=a0+a1*x+a2*(x**2)+a3*(x**3)+a4*(x**4)+a5*(x**5)+a6*(x**6)

ww(j)=at(j)/dsqrt(p)

ot(j)=p+b0+b1*x+b2*(x**2)+b3*(x**3)+b4*(x**4)+b5*(x**5)

ot(j)=ot(j)+b6*(x**6)

at(j)=ww(j)*dcos(ot(j))

ot(j)=ww(j)*dsin(ot(j))

fu(j)=fu(j)+2.d0*at(j)

40 ff(j)=ff(j)+2.d0*ot(j)

fn(j)=dcmplx(fu(j),-ff(j))*h/3.D0

1 fn(j)=0.5d0*ck*fn(j)

 

* Складання матриці h(v)

do 99 i=1,nn

do 99 j=1,nn

l=abs(j-i)+1

fun(i,j)=fn(l)

han1(i,j)=fn(l)

99 jot1(i,j)=fn(l)

 

* Визначення правої частини рівняння F(y)

do 8411 i=1,nn

yyy(i)=yn+t*(i-1)

y=ck*(yyy(i)+t/2.d0)

if(y-3.d0) 842,842,843

842 at(i)=1.+d0*((y/3.d0)**2)+d1*((y/3.d0)**4)+d2*((y/3.d0)**6)

at(i)=at(i)+d3*((y/3.d0)**8)+d4*((y/3.d0)**10)+d5*((y/3.d0)**12)

ot(i)=(2.d0/pi)*dlog(y/2.d0)*at(i)+z0+z1*((y/3.d0)**2)

ot(i)=ot(i)+z2*((y/3.d0)**4)+z3*((y/3.d0)**6)+z4*((y/3.d0)**8)

ot(i)=ot(i)+z5*((y/3.d0)**10)+z6*((y/3.d0)**12)

ot(i)=ot(i)

goto 841

843 at(i)=a0+a1*(3.d0/y)+a2*((3.d0/y)**2)+a3*((3.d0/y)**3)

at(i)=at(i)+a4*((3.d0/y)**4)+a5*((3.d0/y)**5)+a6*((3.d0/y)**6)

at(i)=at(i)/dsqrt(y)

ot(i)=y+b0+b1*(3.d0/y)+b2*((3.d0/y)**2)+b3*((3.d0/y)**3)

ot(i)=ot(i)+b4*((3.d0/y)**4)+b5*((3.d0/y)**5)+b6*((3.d0/y)**6)

ww(i)=at(i)*dcos(ot(i))

ot(i)=at(i)*dsin(ot(i))

at(i)=ww(i)

841 fy(i)=sks*dcmplx(-at(i),ot(i))

rez(i)=fy(i)

8411 u1(i)=fy(i)

 

kl=hpd

if(kl-1)33,33,34

33 an=shr*10.d0+(nn-lp-shr*10.d0)

goto 35

34 an=shr*10.d0+(nn-lp-shr*kl*10.d0)/(kl-1)

 

if(hpd)35,35,98

98 if(((nn-lp-shr*kl*10.d0)/(kl-1))-bda)77,77,35

77 write(1,*)'Неправильно введені дані! Відстань м/д канавками!'

return

 

35 ml=an

an=shr*10.d0

mm=an

 

* Обчислення імпедансу

78 if (hpd)66,66,67

66 do 891 j=1,nn-lp

891 ro(j)=dcmplx(0.d0,ni+0.01*(ll-1))

goto 68

67 do 545 i=1,nn-lp

ro(i)=dcmplx(0.d0,0.d0)

do 87 l=1,kl

do 69 j=1+(l-1)*ml,mm+ml*(l-1)

69 ro(j)=dcmplx(0.d0,ni+0.01*(ll-1))

87 kl=kl

545 m=m

 

* Обчислення значень щільності струму за допомогою рекурентних формул

68 do 8 k=1,nn

dr=ro(k)/(1.d0+ro(k)*jot1(k,k))

do 3 i=1,nn

do 3 j=1,nn

han(i,j)=jot1(i,j)-dr*jot1(i,k)*jot1(k,j)

3 mat(i)=u1(i)-jot1(i,k)*u1(k)*dr

do 11 i=1,nn

ss(i)=mat(i)

do 11 j=1,nn

jot1(i,j)=han(i,j)

11 u1(i)=mat(i)

8 d=d

do 911 i=1,nn

write(5,*)abs(ss(i))

write(6,*)real(ss(i))

911 write(7,*)aimag(ss(i))

* Визначення коефіцієнта придушення й обчислення інтегралів виду I(Y,R),I(Y,0)

vvas=0.D0

vas=0.D0

do 991 i=nn-lp,nn

vvas=vvas+((cdabs(rez(i)))**2)*t

991 vas=vas+((cdabs(ss(i)))**2)*t

vvs=10.d0*dlog10(vvas/vas)

write(*,*)ll,ni+sg*(ll-1),bda

write(4,*)ni+sg*(ll-1)

123 write(1,*)vvs

333 write(3,*)bda

return

end

 

Додаток В

Копія демонстраційного креслення 1

 

 

Рисунок 1-Геометрия задачі

 

де

 

- щільність поверхневого електричного струму;

- функція Ханкеля другого роду нульового порядку з різностним ядром, яке має логарифмічну особливість;

- поверхневий імпеданс;

- кутова частота ();

- абсолютна діелектрична проникність (Ф/м).

 

 

Додаток Г

Копія демонстраційного креслення 2

 

де

- середина –го елементарного інтервалу ;

– число елементарних інтервалів, на які розбито проміжок інтегрування ;

 

де

- нормований сторонній імпеданс на смугах і ;

- комплексна амплітуда щільності електричного струму на смузі d, на якій розташована апертура прийомної антени;

- комплексна амплітуда щільності електричного струму на смузі d при = 0;

 

 

Додаток Д

Копія демонстраційного креслення 3

 

 

Рисунок 2 - Залежність коефіцієнта придушення h від поверхневого імпедансу при ширині

	1 2 3 4 5 6

Рисунок 3 - Залежність модуля щільності струму від при коли імпеданс .

Додаток Е

Копія демонстраційного креслення 4

 

Рисунок 4- Залежність коефіцієнта придушення h від кількості імпедансних смужок шириною (товста лінія), (тонка лінія), (пунктирна лінія), при , коли ширина

	1 2 3 4 5 6 7 8 9 10 11 12 13 14

Рисунок 5 - Залежність модуля щільності струму від , коли кількість імпедансних смужок шириною дорівнює одної (товста лінія), десяти (тонка лінія), при , для ширини

 

Додаток Ж

Копія демонстраційного креслення 5

Таблиця 1 - Значення коефіцієнта придушення η для величіні імпедансу , , , для випадку при

Довжина хвилі ,см	Коефіцієнт придушення
	35,41	-12,29	44,25
	34,51	-10,84	45,97
	33,72	-19,22	48,45
	33,03	-11,25	53,19
	32,42	-24,83	56,57
	31,87	-16,43	-27,04
	31,37	-11,07	-5,71
	30,90	-22,18	-17,53
	30,47	-14,42	-7,39
	30,07	-23,56	-5,97

Таблиця 2 - Значення коефіцієнта придушення η для величіні імпедансу , при , для випадку ,

Довжина хвилі	Коефіцієнт придушення
	11,1245	25,21505
	10,32174	23,01865
	9,676311	20,85377
	9,145582	21,7781
	8,698455	21,85748
	8,313919	21,70971
	7,977633	21,47243
	7,679528	21,19706
	7,412326	20,90591
	7,170615	20,61089