Algoritmi Pentru Rezolvarea Ecuatiilor Integrale

CAPITOLUL I

ECUAȚII INTEGRALE

1.1 INTRODUCERE

Ecuația

(0.1)

în care K(x,s) , f(x) sunt funcții date , este un parametru , iar (x) este o funcție necunoscută , se numește ecuația integrală de speța a doua a lui Fredholm, rezervându-se denumirea de ecuație integrală de speța întâi a lui Fredholm pentru ecuația

(0.2)

Matematicianul suedez Ivar Fredholm (1866-1927) a fost cel dintâi care a făcut o teorie completă a aceste ecuații (1900). În ecuația (0.1) se presupune , în general , că funcția f(x) este continuă pe intervalul [a,b],iar funcția K(x,s) este , de asemenea continuă pe pătratul D definit de inegalitățile

D: a x b a s b

Putând avea discontinuități de-a lungul dreptei x=s.

Probleme de fizică matematică și chiar probleme simple asupra ecuațiilor diferențiale conduc la ecuația integrală (0.1) a lui Fredholm. De exemplu să considerăm ecuația

(0.3)

și să tratăm următoarea problemă la limită asupra ei:

Să se determine soluția ecuației (0.3) care să satisfacă la condițiile bilocale y(a)=0 ,y(b)=0. Se poate proceda în felul următor : se ia ca funcție necunoscută derivata adică

y” =(x) (0.4)

și printr-un calcul ușor se arată că soluția ecuației diferențiale simple (0.4) care satisface la condițiile y(a)=0, y(b)=0 este

(0.5)

unde funcția G(x,s) este definită în modul următor

(0.6)

Într-adevăr soluția generală a ecuației diferențiale (0.4) este

(0.7)

unde C1 și C2 sunt două constante oarecare.

Din prima condiție y(a)=0 rezultă că C2 =0. Din a doua condiție (b)=0, avem

de unde rezultă că

Revenind la ecuația (0.7) , deducem că soluția ecuației diferențiale (0.4) , care satisface la condițiile la limită y(a)=0 ,y(b)=0, este

și aceasta se mai scrie sub forma

adică , avem formula (0.5) unde funcția G(x,s) este dată de formula (0.6).

Scriind că funcția y dată de ecuația (0.5) verifică ecuația (0.3) , suntem conduși la următoarea ecuație pentru funcția necunoscută

(0.8)

unde

K(x,s)=A(x)G(x,s) (0.9)

care este o ecuație integrală a lui Fredholm de forma (0.1) considerată mai sus.

În cazul particular , când în ecuația (0.1) funcția K(x,s) este nulă pentru s>x și continuă pentru a s x b, ecuația se scrie

(0.10)

și se numește ecuația integrală de speța a două a lui Volterra , după numele matematicianului italian Vito Volterra (1860-1940) , care a studiat-o cel dintâi (1896).

Iată de exemplu o problemă care conduce la ecuația integrală (0.10) a lui Volterra:

Să se determine soluția ecuației diferențiale (0.3) care să satisfacă la condițiile y(a)=0 , y’(a)=0 (problema lui Cauchy). Luând ca funcție necunoscută tot derivata y” , așa cum arată ecuația (0.4) , se constată fără greutate că soluția ecuației simple (0.4) , care satisface la condițiile y(a)=0, y’(a)=0, este dată de ecuația

(0.11)

Scriind că funcția (0.11) este soluție a ecuației diferențiale (0.3)suntem conduși la ecuația integrală a lui Volterra (0.10) , în care K(x,s)=A(x)(x-s).

Deoarece ecuația (0.10) a lui Volterra este un caz particular al ecuației integrale (0.1) a lui Fredholm , în acest capitol vom studia numai ecuația integrală a lui Fredholm de speța a doua.

Printre fondatorii teoriei ecuației integrale liniare , în afară de Volterra și Fredholm , mai cităm pe David Hilbert(1862-1943) și Erhart Schmitdt (1876-1958). Mai cităm pe matematicianul român Traian Lalescu , care în teza sa de doctorat “ Sur l’equation de Volterra”, susținută la Paris (1908) , a aplicat cel dintâi metoda aproximațiilor succesive la integrarea ecuației lui Volterra. Tot Lalescu a scris cea dintâi carte din lume despre ecuațiile integrale care a apărut la București în limba română (1911), și un an mai târziu la Paris în limba franceză (1912). Această carte a jucat un rol important în dezvoltarea teoriei ecuațiilor integrale, carte care și azi este citită cu mare interes.

1.2 APLICAREA METODEI APROXIMAȚIILOR SUCCESIVE

LA INTEGRAREA ECUAȚIEI INTEGRALE A LUI FREDHOLM

Să considerăm ecuația integrală de speța a doua a lui Fredholm

(1.1)

și să formulăm condițiile în care se studiază această ecuație. Vom face următoarele ipoteze.

1.Funcția K(x,s) numită nucleul ecuației integrale este o funcție complexă continuă pe domeniul D.

2.Funcția f(x) numită termenul liber al ecuației integrale este o funcție complexă , continuă pe intervalul [a,b].

3. este un parametru, în general număr complex.

Se caută soluția ecuației integrale (1.1) în mulțimea funcțiilor continue pe intervalul [a,b] sau în spațiul funcțiilor continue pe intervalul [a,b] , care se notează cu C[a,b]sau mai scurt cu C.

Teorema de existență și unicitate. În condițiile formulate în (1.1) și , unde m este o margine superioară a modulului nucleului K(x,s) pe domeniul D, ecuația integrală (1.1) are o soluție continuă pe intervalul [a,b], unică. Ea poate fi obținută cu metoda aproximațiilor succesive.

Pentru a demonstra această teoremă , să observă că la ecuația (1.1) și în condițiile specificate mai sus se poate aplica principiul contracției.

Să introducem operatorul

(1.2)

unde

(1.3)

este un operator numit operatorul lui Fredholm.

Ținând seama de aceasta , ecuația integrală (1.1) se scrie sub forma

A()= (1.4)

Deoarece funcția f(x) este continuă pe intervalul [a,b], iar operatorul lui Fredholm K face să corespundă la orice funcție continuă pe intervalul [a,b], operatorul A() transformă elementele spațiului C în elemente ale aceluiași spațiu. Distanța dintre elementele și ale spațiului C se definește prin

(1.5)

Observăm că

și prin urmare avem

Însă și introducând numărul

(1.6)

vom avea PM(b-a), unde m este o margine superioară a lui pe domeniul D. Rezultă că

și prin urmare

(1.7)

unde

(1.8)

Atunci , aplicând principiul contracției , ecuația (1.4) are o soluție unică și ea poate fi obținută cu metoda aproximațiilor succesive. Convergența șirului (n) către trebuie înțeleasă în sensul metricei spațiului C , introdusă rin formula (1.5) adică este convergența uniformă.

Puterile succesive ale operatorului K() și nucleii iterați. Din aplicarea principiului contracției rezultă că șirul

converge uniform către soluția a ecuației integrale (1.1) când . De aici rezultă că scrie

(1.9)

Introducând funcțiile

)

vom putea scrie

adică

de unde rezultă că

Seria (1.9) devine:

cu

(1.10)

prin urmare soluția (x) a ecuației integrale (1.1) se dezvoltă într-o serie după puterile lui ,

(1.11)

unde coeficienții sunt funcții continue de x pe intervalul [a,b], date de formulele (1.10) și care converg absolut și uniform pe intervalul [a,b] când .

Definiție. Puterea a n-a a operatorului Kf , care se notează cu Knf este dată de formula

Knf=K(Kn-1f), (1.12)

pentru n2 și K1f=Kf.

Definiție. Nucleul iterat de ordinul n , care se notează cu Kn(x,s), este dat de formula

(1.13)

pentru n2 și K1(x,s)=K(x,s).

Pentru n=2 formulele (1.12) și (1.13) ne arată că avem

(1.14)

iar pentru n>2 vom avea

(1.15)

Observăm din formula (1.13) că nucleul Kn(x,s) se mai scrie

(1.16)

și din această formulă mai rezultă că putem scrie formula (1.13) sub forma

(1.17)

Operatorii Knf astfel introduși au proprietatea

(1.18)

care rezultă din proprietatea nucleilor iterați

(1.19)

Aceasta se stabilește în modul următor. Din formula (1.16) se deduce că avem

Deoarece prima paranteză este Kn(x,t) , iar a doua este Km(x,t), suntem conduși la formula (1.19).

Pentru a demonstra formula (1.18) , ne servim de formula (1.15) și vom avea

adică

de unde rezultă formulele (1.18).

Cu introducerea nucleilor iterați , formulele (1.10) se pot transforma și scrie sub forma

(1.20)

pentru n=1,2,…….

Pentru n=1, formula (1.20) este evidentă: ea rezultă din primele două formule (1.10). În general formula ( 1.20) se demonstrează cu metoda inducției complete. Să admitem că pentru indicele n avem

(1.21)

ceea ce este adevărat pentru n=1. Atunci , conform formulei (1.10) pentru indicele n+1 și formulei (1.21) avem

și aceasta se mai scrie sub forma

conform formulei (1.13) pentru indicele n+1. Astfel formulele (1.20) sunt demonstrate.

Forma soluției ecuației integrale (1.1). S-a arătat în (1.3)că, în catul când soluția ecuației integrale a lui Fredholm (1.1) se dezvoltă într-o serie (1.11) după puterile lui . Să notăm cu n(x) suma primilor n+1 termeni ai acestei serii , adică

(1.22)

Ținând seama de formulele (1.20) , aceasta se mai scrie

n(x,s;))f(s)ds (1.23)

unde s-a notat

n(x,s;)=K(x,s)+ K2(x,s)+……. n-1Kn(x,s) (1.24)

Observăm că suma n(x,s;) este suma primilor n termeni ai seriei

K(x,s)+ K2(x,s)+……. n-1Kn(x,s)+….. (1.25)

Să demonstrăm că această serie este absolut uniform convergentă pe pătratul D și că suma ei este soluția ecuației integrale

(1.26)

pe care o mai putem scrie sub forma

(1.27)

unde s-a notat cu A[(x,s)] operatorul definit de membrul al doilea al ecuației (1.26).

La ecuația (1.26) sau(1.27) se poate aplica principiul contracției în același mod cum s-a procedat și la ecuația integrală (1.4) , demonstrându-se teorema de existență și unicitate. Este suficient pentru aceasta să se definească distanța (,) dintre două funcții (x,s) și (x,s) prin formula

Dacă soluția ecuației (1.26) există și este unică. Deoarece soluția se poate obține cu metoda aproximațiilor succesive , rezultă, refăcând raționamentele care ne-au condus la serie (1.11) , că ea poate dezvolta în serie după puterile lui , absolut și uniform convergentă pe pătratul D,

(1.28)

unde avem, analog cu formulele (1.10)

(1.29)

Comparând formulele (1.29) cu formulele (1.13) , deducem că

(1.30)

astfel încât seria (1.28) coincide cu seria (1.25) și , prin urmare am demonstrat că seria (1.25) converge absolut uniform pe pătratul D , când

Se notează de obicei

(x,s;)=K(x,s)+ K2(x,s)+…….+ n-1Kn(x,s)+… (1.31)

și funcția (x,s;) se numește nucleul rezolvant al ecuației integrale (1.1) a lui Fredholm.

Suma n(x) tinde uniform către soluția (x) a ecuației integrale (1.1) când n, iar suma (x,s;) tinde de asemenea către nucleul rezolvant (x,s;) când n.

Rezultă atunci că dacă în formula (1.23) facem pe n și observăm că

n(x,s;)f(s)ds(x,s;)f(s)ds ,

deducem că

(x,s;)f(s)ds. (1.32)

Această formulă rezolvă ecuația lui Fredholm în cazul când și se vede bine rolul nucleului rezolvant.

Identitățile pe care le verifică nucleul rezolvant. S-a văzut în alineatul precedent că nucleul rezolvant (x,s;) este soluția ecuației integrale (1.26) ceea ce înseamnă că avem prima identitate

(x,s;) =K(x,s) +K(x,t) (x,s;)dt, (1.33)

pe care o verifică nucleul rezolvant.

Să transformăm funcția n(x,s;) dată de formula (1.24) cu ajutorul formulelor (1.17). Schimbând pe n cu n+1 vom avea

n+1(x,s;)=K(x,s)+ +

+

care se mai scrie sub forma

n+1(x,s;)=K(x,s)+

adică avem

n+1(x,s;)=K(x,s)+ n(x,s;)K(t,s)dt rezultă, refăcând raționamentele care ne-au condus la serie (1.11) , că ea poate dezvolta în serie după puterile lui , absolut și uniform convergentă pe pătratul D,

(1.28)

unde avem, analog cu formulele (1.10)

(1.29)

Comparând formulele (1.29) cu formulele (1.13) , deducem că

(1.30)

astfel încât seria (1.28) coincide cu seria (1.25) și , prin urmare am demonstrat că seria (1.25) converge absolut uniform pe pătratul D , când

Se notează de obicei

(x,s;)=K(x,s)+ K2(x,s)+…….+ n-1Kn(x,s)+… (1.31)

și funcția (x,s;) se numește nucleul rezolvant al ecuației integrale (1.1) a lui Fredholm.

Suma n(x) tinde uniform către soluția (x) a ecuației integrale (1.1) când n, iar suma (x,s;) tinde de asemenea către nucleul rezolvant (x,s;) când n.

Rezultă atunci că dacă în formula (1.23) facem pe n și observăm că

n(x,s;)f(s)ds(x,s;)f(s)ds ,

deducem că

(x,s;)f(s)ds. (1.32)

Această formulă rezolvă ecuația lui Fredholm în cazul când și se vede bine rolul nucleului rezolvant.

Identitățile pe care le verifică nucleul rezolvant. S-a văzut în alineatul precedent că nucleul rezolvant (x,s;) este soluția ecuației integrale (1.26) ceea ce înseamnă că avem prima identitate

(x,s;) =K(x,s) +K(x,t) (x,s;)dt, (1.33)

pe care o verifică nucleul rezolvant.

Să transformăm funcția n(x,s;) dată de formula (1.24) cu ajutorul formulelor (1.17). Schimbând pe n cu n+1 vom avea

n+1(x,s;)=K(x,s)+ +

+

care se mai scrie sub forma

n+1(x,s;)=K(x,s)+

adică avem

n+1(x,s;)=K(x,s)+ n(x,s;)K(t,s)dt (1.34)

În formulele (1.34) să facem pe n și să observăm că membrul întâi n+1(x,s;) tinde către nucleul rezolvant (x,s;), iar în membrul al doilea

n(x,t;)K(t,s)dt(x,t;)K(t,s)dt

Deducem astfel o a doua identitate

n+1(x,s;)=K(x,s)+ n(x,t;)K(t,s)dt (1.35)

pe care o verifică nucleul rezolvant.

Ecuația integrală asociată. Ecuația integrală

(1.36)

unde nucleul este conjugatul nucleului care se obține din k(x,s) permutând pe x cu s, iar este conjugatul lui se numește ecuația integrală a ecuației (1.1). Nucleul

(1.37)

este nucleul asociat nucleului lui K(x,s).

1.Între nucleii iterați și avem aceeași relație ca între nucleii K*(x,s) și K(x,s), adică

(1.38)

Într-adevăr , aplicând formulele (1.16) , avem

și , prin urmare , formulele (1.38) sunt dovedite.

2.Din formulele (1.38) rezultă că nucleul rezolvant al ecuației (1.36) (x,s,), este dat de formula

(x,s,)= (1.39)

Într-adevăr, avem

(x,s,)=

și, ținând seama de formula (1.38) , obținem formula (1.39).

3.Soluția ecuației (1.36) pentru este dată de formula (1.32) și ținând seama de formula (1.39), vom avea

(1.40)

Am obținut, prin urmare următoarea

Teoremă. Dacă , ecuația integrală (1.1) a lui Fredholm are o soluție unică dată de formula (1.32) , iar ecuația integrală asociată (1.36) are , de asemenea, o soluție unică dată de formula (1.40).

Prin urmare, rezolvarea ecuației integrale a lui Fredholm este făcută când Însă , problema rezolvării ecuației lui Fredholm se pune și pentru

1.3 ECUAȚIA INTEGRALĂ A LUI VOLTERRA

Să considerăm ecuația integrală a lui Fredholm

(2.1)

și să presupunem că nucleul K(x,s) este o funcție continuă în domeniul asxb, că este , de asemenea , o funcție continuă și mărginită în domeniul axsb , dar este discontinuă de-a lungul dreptei x=s.

Se arată cu ușurință că se poate aplica principiul contracției la ecuația (2.1).

Rezultă atunci că ecuația integrală (2.1) are o soluție unică , pe care o putem obține cu metoda aproximațiilor succesive , dacă unde pe pătratul D.

Ecuația integrală a lui Volterra , definită la începutul acestui capitol adică ecuația

(2.2)

unde K(x,s) este un nucleu continuu pe triunghiul T definit de inegalitățile asxb, intră în clasă de ecuații integrale (2.1).

Ea are deci o soluție unică , pe care o putem obține cu metoda aproximațiilor succesive. Însă, la ecuația lui Volterra avem următorul rezultat important , dat de următoarea:

Teoremă. Soluția ecuației integrale (2.2) a lui Volterra se poate dezvolta într-o serie întreagă după puterile lui , uniform convergentă pe intervalul [a,b] și oricare ar fi .

Într-adevăr , din teorema de existență și unicitate rezultă că soluția ecuației (2.2) a lui Volterra se dezvoltă într-o serie întreagă după puterile lui

(2.3)

unde

Notând

(2.4)

se arată cu ușurință că avem

(2.5)

de unde rezultă convergența absolută și uniformă a seriei (2.3) oricare ar fi . Soluția (x) este dată de formula lui Volterra

(x,s;)f(s)ds, (2.6)

unde (x,s;) este nucleul retolvant , definit de seria

(x,s;)= (2.7)

unde Kn(x,s) sunt nucleii iterați

(2.8)

Seria (2.7) este absolut uniform convergentă pe triunghiul T , oricare ar fi . Prin urmare , la ecuația integrală a lui Volterra nu poate fi vorba de valorii proprii și de funcții proprii.

CAPITOL II

2.1 METODE DE REZOLVARE A

ECUAȚIILOR INTEGRALE

INTRODUCERE

Mulți oameni , altminteri destul de puțini din punct de vedere numeric, își imaginează că soluția numerică a ecuațiilor integrale are o argumentație încâlcită deoarece, până de curând, această soluție nu a fost explicată în cărțile școlare / manuale de analiză numerică. De fapt, există o preocupare în continuă dezvoltare despre soluția numerică a ecuațiilor integrale , câteva monografii au apărut deja.

Unul dintre motivele slabei activități în acest domeniu ar fi acela că există multe tipuri de ecuații, fiecare având diverse capcane posibile ; adesea fiind propuși mai mulți algoritmi pentru un singur caz.

Există o strânsă corespondență între ecuațiile liniare integrale , care definesc relații liniare și integrale între funcții definite pe un spațiu de dimensiune infinită și vechile ecuații simple și liniare , care denumesc relații asemănătoare între vectori într-un spațiu vectorial finit-dimensional. Deoarece această corespondență stă la baza multor algoritmi computaționali merită să recapitulăm criteriile de calificare a ecuațiilor integrale.

Ecuațiile Fredholm implică integrale definite cu limite superioare și inferioare. O ecuație Fredholm neomogenă de primul tip are forma:

(0.1)

Aici f(t) e funcția necunoscută care trebuie rezolvată , în timp ce g(t) e cunoscută ca "partea din dreapta". (În ecuațiile integrale pentru niște motive foarte ciudate familiara "partea dreaptă" e scrisă convențional în stânga). Funcția de două variabile R(t,s) se numește nucleul. Ecuația (0.1) e analoagă cu ecuația matrice:

(0.2)

a cărei soluție e f=K-1g, unde R-1 e matricea inversă. Ca și ecuația (0.2) ecuația (0.1) are soluție unică ori de câte ori g 0 (cazul omogen cu g=0 nu e niciodată de ajutor) și K e inversabil. Oricum așa cum vom vedea , această ultimă condiție e mai des o excepție decât o regulă.

Analogul problemei cu valoarea proprie finit-dimensională

(K-1)f=g (0.3)

e numită o ecuație Fredholm de tipul 2 , scrisă de obicei

(0.4)

de asemenea, convențiile de notare nu corespund în mod exact: în ecuația (0.4) este 1/ în (0.3), în timp ce g este -g/. Dacă g=0 atunci ecuația e numită omogenă. Dacă nucleul K(t,s) e mărginit , atunci , ca și ecuația (0.3), ecuația (0.4) are proprietatea că forma sa omogenă are soluții pentru o mulțime infinită , un interval de dimensiune cel mult numărabilă =n , n=1,2,3…, valori proprii. Soluțiile corespunzătoare fn(t) sunt funcții cu valori proprii. Valorile proprii sunt reale dacă nucleul e simetric.

În cazul neomogen când g#0 ecuațiile (0.3) și (0.4) sunt rezolvabile cu excepția cazului când (sau ) este o valoare proprie deoarece atunci operatorul integral (sau matriceal) este singular. În ecuațiile integrale această diferență este numită "alternativa Fredholm".

Ecuațiile Fredholm de tipul întâi sunt deseori foarte prost definite.

Aplicarea nucleului la o funcție este în general o operație omogenă fără întreruperi , astfel încât soluția care cere inversarea operatorului , va fi foarte sensibilă la micile schimbări sau erori ale variabilelor. De fapt omogenizarea pierde des informații și nu există nici o cale de a le recupera printr-o operație inversă. S-au dezvoltat metode specializate pentru astfel de ecuații , ecuații denumite deseori "problemele inversei".

În general , o metodă trebuie să adauge informației niște cunoștințe anterioare asupra naturii soluției. Apoi , această informație anterioară e folosită într-un fel sau altul , pentru a recupera informația pierdută. Vom introduce aceste tehnici în capitolul IV.

Ecuațiile neomogene Fredholm de tipul al doilea sunt prost definite mult mai rar. Ecuația (0.4) poate fi scrisă astfel

(0.5)

unde (t-s) e o funcție delta Dirac )și unde am schimbat de la la analogul pentru claritate. Dacă este destul de mare ca dimensiune , atunci ecuația (0.5) este , de fapt, diagonal dominantă și deci bine definită. Doar dacă este mic atunci ne întoarcem la cazul de proastă definire.

Ecuațiile omogene Fredholm de tipul al doilea sunt de asemenea prost poziționate dar nu în mod special. Dacă R este un operator omogen atunci va înlocui mulți f ' s care tind spre zero sau aproape de zero ; va exista astfel un mare un mare număr de valori proprii degenerate sau aproape degenerate în jurul =0() , dar acesta nu va cauza nici o dificultate de evaluare. De fapt putem vedea acum că mărimea lui care este necesară pentru a putea salca ecuația neomogenă (0.5) de o proastă definire mult mai puțin decât cea cerută de metoda dominării diagonală.

O dată se termenul schimbă toate valorile proprii , este îndeajuns încât să fie destul de mare pentru a schimba unicitatea intervalului operatorului de unicitate al valorilor proprii apropriate și depărtate de zero , astfel încât operatorul care rezultă de aici devine inversabil (cu excepția valorilor proprii discrete).

Ecuațiile Volterra sunt un caz special al ecuațiilor Fredholm cu K(t,s)=0 pentru s>t. Dând la o parte partea inutilă a integrării , ecuațiile Volterra sunt scrise într-o formă în care limita superioară a integrării este variabila independentă t.

Ecuația Volterra de primul tip

(0.6)

are ca analog ecuația matrice (acum scrisă în componente)

(0.7)

Prin compararea cu ecuația (0.2) vedem că ecuațiile Volterra corespund unei matrici K care este inferior (stânga) triunghiulară cu toate valorile de deasupra diagonalei egale cu zero. Se știe că astfel de ecuații matrice sunt rezolvabile în mod obișnuit printr-o substituire anterioară.

Tehnicile pentru rezolvarea ecuațiilor Volterra sunt de asemenea directe. Atunci când "zgomotul" de măsurare experimental nu este dominant, ecuațiile Volterra de primul tip tind să nu fie prost definite ; limita superioară a integralei introduce un pas clar care strică orice proprietăți de omogenizare ale nucleului.

Ecuația Volterra de tipul doi este scrisă:

(0.8)

a cărei matrice analoagă este ecuația

(K-1)f=g (0.9)

cu k inferior triunghiulară. Motivul pentru care nu există în aceste ecuații este acela că (i) în cazul neomogen diferit de zero acesta poate fi inclus în K , în timp ce (ii) în cazul omogen (g=0) este o teoremă că ecuația de tipul doi cu nuclee mărginite nu au valori omogene cu funcții omogene dublu integrabile.

Ne-am specializat definițiile în cazul ecuațiilor liniare integrabile.

Componenta într-o versiune neliniară a ecuațiilor (0.1) sau (0,6) ar deveni K(t,s,f(s)) în loc de K(t,s)f(s); o versiune neliniară a ecuațiilor (0.4) sau (0.8) ar avea o componentă K(t,s,f(t),f(s)). Ecuațiile neliniare Fredholm sunt considerabil mai complicate decât echivalentele lor liniare. Din fericire ele nu apar atât de des în practică și în general vorbind le vom ignora în acest capitol. Din contră , rezolvarea ecuațiilor neliniare Volterra implică de obicei doar o mică modificare a algoritmului pentru ecuații liniare , așa cum vom vedea în continuare.

Aproape toate metodele pentru rezolvarea ecuațiilor integrale numerice se folosesc de regulă cvadraturi , de regulă cvadraturile Gaussian. În secțiunile următoare vom discuta mai întâi ecuațiile Fredholm de tipul doi cu nuclee omogene(cap. 1). Regulile de cvadratură specială sunt și ele discutate , dar ne vom ocupa de sistemele de ecuații bine definite. Apoi ne întoarcem la ecuațiile Volterra(cap.2) și aflăm că metodele simple și directe sunt în general satisfăcătoare pentru aceste ecuații.

În acest capitol discutăm despre cum să procedăm în cazul nucleelor singulare îndreptându-ne atenția în mod special mai pe larg spre ecuațiile Fredholm (de tipurile 1 și 2).

Cazurile unice cer reguli de cvadratură speciale , dar uneori sunt partea bună a unei nenorociri deoarece poate strica omogenitatea unui nucleu și să creeze probleme a definirii nucleului.

Unele subiecte ca ecuații integrale diferențiale trebuiesc declarate în afara scopului nostru.

Se înțelege că această lucrare abia abordează pe larg câteva din cele mai elementare metode incluse în acest subiect complicat.

2.2. ECUAȚIILE FREDHOLM DE TIPUL DOI

Se cere o soluție numerică pentru f(t) în ecuația:

(1.1)

Metoda descrisă , una foarte elementară , este numită metoda Nystrom și cere alegerea unei reguli de cvadratură pentru aproximarea soluției.

(1.2)

Aici mulțimea {wj} are ca elemente valorile cvadraturii în timp ce punctele {sj} sunt abscisele.

Ce cvadratură să folosim? E posibil să rezolvăm ecuațiile integrale cu sisteme inferioare cum ar fi repetarea metodei trapezului sau regulile lui Simpson. Oricum vom vedea că găsirea metodei pentru aflarea soluției implică operații O(N3) și astfel cele mai eficiente metode tind să folosească sisteme superioare de cvadratură pentru a-l menține pe N cât mai mic. Pentru probleme omogene și nesingulare , nimic nu e mai bun decât cvadratura Gaussiană. Pentru nuclee neomogene sau singulare , a se vedea capitolul 3.

Delves și Mohamed[1] au investigat metode mai complicate decât metoda Nystrom. Pentru ecuații Fredholm directe de tipul doi au ajuns la concluzia "…câștigătoarea acestui concurs a fost procedeul Nystrom … cu regula celor N puncte Gauss-Legendre. Acest procedeu este foarte simplu. … Aceste rezultate sunt suficiente pentru a face analiza numerică 'să plângă' ".

Dacă aplicăm cvadratura (1.2) la ecuația (1.1) avem

(1.3)

Evaluăm ecuația (1.3) la punctele de cvadratură

(1.4)

Fie fi vectorul f(ti) , gi vectorul g(ti) , Kij matricea K(ti,sj) și definim (1.5)

Apoi în matricea ecuației de matrici (1.4)

(1.6)

Aceasta e o mulțime N de ecuații algebrice liniare cu N necunoscute care poate fi rezolvată prin tehnicile standard de descompunere triunghiulară – aici este locul unde numărul de operații de O(N3) își face apariția. Soluția e de obicei bine definită , doar dacă e foarte aproape de valorile proprii.

Dacă s-a găsit soluția la punctele de cvadratură {ti}, cum găsim soluția la orice alt punct. În nici un caz nu trebuie folosită interpolarea polinomială , deoarece aceasta distruge toată acuratețea pentru care s-a muncit din greu s-o obținem. Observația cheie a lui Nystrom a fost aceea că trebuie folosită ecuația (1.3) ca pe o formulă de interpolare , menținând precizia soluției. Vom da aici două procedee pentru ecuațiile liniare Fredholm de tipul al doilea. Procedeul fred 2 pregătește ecuația (1.6) și apoi o rezolvă prin descompunerea LU care face apel la procedurile ludcmp și lubksb. Cvadratura Gauss-Legendre este implementată prin obținerea mai întâi absciselor cu o apelare spre gauleg. Procedura fred 2 cere să se asigure o funcție externă care să returneze g(t) și alta care să returneze Kij. Apoi se returnează soluția f la punctele de cvadratură. De asemenea se returnează punctele și valorile de quadratură. Acestea sunt folosite de a doua procedură fredin pentru a realiza interpolarea Nystrom a ecuației (1.3) și returnează valoarea lui f la orice punct din intervalul {a,b}.

#include "nrutil.h"

void fred 2 (int n, float a, float b, float t[], float f[] , float w[] ,

float (*g)(float) , float(*ak)(float,float))

Rezolvă o ecuație liniară Fredholm de tipul al doilea. Ca date de intrare a,b sunt limitele integrării și n este numărul de puncte folosite în cvadratura Gaussiană g și ak , acestea fiind funcții externe furnizate de utilizator , care returnează g(t) și K(t,s). procedura returnează vectorii t[1..n] și f[1..n] conținând abscisele ti ale cvadraturii Gaussiene și soluția f la aceste abscise. De asemenea este returnat și vectorul w[1..n] al valorilor Gaussiene care se folosesc cu procedura de interpolare Nystrom fredin

{

void gauleg(float x1, float x2, float x[], float w[], int n);

void lubksb(float **a, int n, int *indx, float b[]);

void ludcmp(float **a, int n, int *indx, float *d);

int i,j, *indx;

floar d, **amk;

indx=ivector(1,n);

omk=matrix(1,n,1,n);

gauleg(a,b,t,w,n); \se înlocuiește cu altă procedură dacă nu folosim

for (i=1;i<=n;i++); { \cvadratura G-L , formăm 1-

for (j=1;j<=n; j++)

ommk[i][j]=float (i == j)-(*ak)(t[i],t[j])*w[j];

f[i]=(*g)(t[i]);

}

ludcmp(omk,n,indx,&d); \rezolvă ecuația liniară

lubksb(omk,n,indx,f);

free_matrix(omk,1,n,1,n);

free_ivector(indx,1,n);

}

floatfredin(float x,int n,float a,floatb,float t[], float f[],

float w[], float (*g)(float) ,float(*ak)(float,float))

Fiind dați vectorii t[1..n] și w[1..n] conținând abscisele și valorile cvadraturii Gaussiene și dându-se vectorul soluție f[1..n] din fred 2, această funcție returnează valoarea lui f în punctul x folosind formulade interpolare Nystrom. Ca date de intrare a și b sunt limitele integrării și n e numărul de puncte folosite în cvadratura Gaussiană. g și ak sunt funcții externe furnizate de utilizator care returnează g(t) și respectiv K(t,s).

{

int i;

float sum=0.0;

for (i=1;i<=n;i++) sum += (*ak)(x,t[i])*w[i]*f[i];

return (*g) (x)+sum;

}

Un dezavantaj al metodei bazate pe cvadratura Gaussiană este acela că nu există o cale simplă pentru estimarea erorii din rezultat. Cea mai practică metodă este să creștem N cu 50% , să zicem , tratăm diferența dintre cele două estimări ca fiind eroarea din rezultatul obținut cu valoarea mai mare a lui N.

Să ne întoarcem acum la soluțiile ecuațiilor omogene. Dacă dăm valoarea =1 și g=0 , atunci ecuația (1.6) devine o ecuație de valoare ? eigen standard

(1.7)

pe care o putem rezolva cu orice procedeu convenabil de matrice cu valori proprii. De notat că dacă problema noastră inițială are un nucleu simetric , atunci matricea K este simetrică. Oricum , odată ce valorile wj nu sunt egale pentru majoritatea cvadraturilor , matricea (ec. 1.5) nu este simetrică. Problema matricei de valori proprii este mai ușoară pentru matricile simetrice și astfel vom simetriza matricea dacă e posibil.

Dacă valorile se dovedesc a fi pozitive (așa cum sunt pentru cvadratura Gaussiană ) putem defini matricea diagonală D=diag(wj) și rădăcina sa pătrată D1/2=diag(). Atunci ecuația (1.7) devine

KDf=f

Înmulțind cu D1/2 obținem

(D1/2KD1/2)h=h (1.8)

unde h=D1/2f. ecuația (1.8) are acum forma unei probleme de valoare proprie simetrică.

Soluția ecuațiilor (1.7) și (1.8) va da în general N valori proprii , unde N este numărul punctelor de cvadratură folosite. Pentru nucleele dublu integrabile acestea vor da o aproximare bună celor mai mici N valori proprii ale ecuațiilor integrale. Nucleele de rang finit (numite de asemenea nuclee degenerate sau separabile) au doar un număr finit de valori proprii diferite de zero (posibil nici una). Se poate diagnostica această situație printr-un fascicul de valori proprii care sunt zero la precizia calculatorului. Numărul de valori proprii nenule va rămâne constant pe măsură ce N crește pentru a le îmbunătăți precizia. Este nevoie de puțină atenție la: un nucleu nedegenerat poate avea un număr infinit de valori proprii care au un punct de acumulare la =0. Se deosebesc cele două cazuri prin comportamentul soluției pe măsură ce N crește. Dacă suspectăm un nucleu degenerat , vom fi capabili să rezolvăm problema prin tehnicile analitice descrise în toate manualele.

2.3 ECUAȚII VOLTERRA

Să ne întoarcem la ecuațiile Volterra al căror prototip este ecuația Volterra de tipul doi:

(2.1)

Majoritatea algoritmilor ecuațiilor Volterra merg de la t=a construind soluția pe parcurs. Din acest punct de vedere se aseamănă nu numai cu substituirea anterioară (așa cum am mai discutat în Introducere), dar și cu problemele de valoare inițială pentru ecuațiile diferențiale obișnuite. De fapt , mulți algoritmi pentru ODEs au corespondenți pentru ecuațiile Volterra. Cea mai ușoară cale de procedare este rezovarea ecuației cu o (celulă,plasă,mesh) cu un interval uniform:

ti=a+ih, i=0,1,…,N, (2.2)

Pentru aceasta trebuie aleasă o cvadratură. Pentru o celulă uniformă , cea mai simplă schemă este metoda trapezelor:

(2.3)

Astfel metoda trapezelor pentru ecuația (2.1) este:

, i=1,…,N (2.4)

(Pentru ecuațiile Volterra de tipul 1 , primul 1 din stânga va fi absent și g va avea semn opus , cu schimbările directe corespunzătoare în restul discuției).

Ecuația (2.4) e o formulă explicită care dă soluția în O(N2) operații.

Spre deosebire de ecuațiile Fredholm nu e necesară rezolvarea unui sistem de ecuații liniare. Ecuațiile Volterra implică de obicei mai puțină muncă decât ecuațiile Fredholm corespunzătoare care , așa cum am văzut, implică inversarea sistemelor liniare , care câteodată pot fi chiar mari.

Eficiența rezolvării ecuației Volterra e într-un fel contrabalansată de faptul că sistemele de ecuații apar mult mai frecvent în practică. Dacă interpretăm ecuația (2.1) ca o ecuație vector pentru m vectori de funcție f(t) atunci nucleul R(t,s) este o matrice mm. Ecuația (2.4) trebuie înțeleasă ca o ecuație vectorială. Pentru fiecare i trebuie să rezolvăm mulțimea mm de ecuații liniare algebrice prin eliminarea Gaussiană.

Procedura voltra de mai jos implementează acest algoritm. Trebuie să dăm o funcție externă care returnează k funcții ale vectorului g(t) în punctul t, și alta care returnează elementul (k,l) al matricei R(t,s) în (t,s). Procedura voltra returnează atunci vectorul f(t) la distanțe egale ti.

#include “nrutil.h”

void voltra (int n,int m, float t0, float h, float *t, float **f,

float (*g)(int,float),float(*ak) (int,int,float,float))

Rezolvă o mulțime de m ecuații Volterra liniare de tipul 2 folosind regula trapezului extinsă. La intrare t0 este punctul de pornire al integrării și n-1 este numărul de pași ai mărimii h de luat. g(t) e o funcție externă dată de utilizator care returnează gk(t) , în timp ce ak(k,l,t,s) e o altă funcție externă dată de utilizator care returnează elementele (k,l) al matricei K(t,s). Soluția este returnată în f[1..m][1..n] cu absicesle corespunzătoare în t[1..n].

{

void lubksb(float **a, int n, int *indx, float b[]);

void ludcmp(float **a, int n, int *indx, float *d);

int i,j,k,l ,*indx;

float d,sum,**a,*b;

indx=ivector(1,m);

a=matrix(1,m,1,m);

b=vector(1,m);

t[1]=t0;

for (k=1;k<=m;k++) f[k][1]=(*g)(k,t[1]);

for (i=2; i<=m ; i++) {

t[i]=t[i-1]+h;

for (k=1;k<=m;k++) {

sum=(*g)(k,t[i]);

for (l=1;l<=m;l++) {

sum+=0.5*h*(*ak)(k,l,t[i], t[1])*f[1][1];

for (j=2; j<i; j++)

sum+=h*(*ak)(k,l,t[i], t[j])*f[1][1];

a[k][l]=(k==1)-0.5*h*(*ak)(k,l,t[i], t[i]);

}

b[k]=sum;

}

ludcmp(a,m,indx,&d);

lubksb(a,m,indx,b);

for (k=1;k<=m;k++ f[k][i]=b[k];

}

free_vector(b,1,m);

free_matrix(a,1,m,1,m);

free_ivector(indx,1,m);

}

Pentru ecuația Volterra neliniară, ecuația (2.4) acceptă produsul Kii(fi) și de asemenea pentru alte două produse ale lui K și ale lui f. Astfel

pentru fiecare i rezolvăm o ecuație liniară pentru fi cu partea dreaptă cunoscută. Metoda lui Newton cu o presupunere inițială a lui fi-1 de obicei merge foarte bine cu condiția ca mărimea pasului să nu fie prea mare.

Metodele de ordin superior pentru rezolvarea ecuațiilor Volterra , sunt în opinia noastră , nu prea importante ca pentru ecuațiile Fredholm , căci ecuațiile Volterra sunt relativ mai ușor de rezolvat. Oricum , există o literatură în continuă creștere pe acest subiect. Dar apar câteva dificultăți. La început orice metodă care atinge un grad superior prin operarea asupra câtorva puncte de cvadratură , simultan va necesita o metodă specială pentru a începe , când valorile de la primele câteva puncte nu sunt cunoscute încă.

În al doilea rând , regulile stabilite de cvadratură pot ridica instabilități neașteptate în ecuații integrale. De exemplu , să presupunem că încercăm să înlocuim regula trapezului în algoritmul de deasupra cu regula lui Simpson. Regula lui Simpson de obicei integrează peste un interval 2h , așa că putem obține ușor valorile funcției chiar la punctele celulei. Pentru punctele de celulă impare , putem încerca să anexăm un tabel al regulii trapezului. Dar la care final al integrării al trebui să-l anexăm? Putem face un pas al regulii trapezului urmat de toată regula lui Simpson , sau regula lui Simpson cu un pas al regulii trapezului la sfârșit. În mod surprinzător , prima schemă este instabilă, în timp ce adoua este bună.

O simplă apropiere care poate fi folosită cu metoda trapezului dată mai sus este extrapolarea Richardson: calculează soluția cu mărimea pasului h și h/2. Apoi presupunând că eroarea este proporțională cu h2 calculează:

FE= (2.5)

Această procedură poate fi repetată ca și integrarea Romberg. Consensul general este acela că cea mai bună metodă din cele de ordin superior este metoda block-by-block (oprire după oprire).

Alt subiect important este folosirea metodelor cu mărimea pasului variabilă care sunt mult mai eficiente dacă sunt caracteristici exacte în K sau f. metodele cu mărimea pasului variabilă sunt ceva mai complicate decât omoloagele lor pentru ecuații diferențiale.

ECUAȚII INTEGRALE CU NUCLEE

SINGULARE UNICE

Multe ecuații integrabile au particularități ori în nucleu ori în soluție ori în amândouă. O simplă metodă de cvadratură va arăta o convergență necorespunzătoare cu N dacă aceste particularități sunt ignorate. Există uneori un truc prin care particularitățile sunt cel mai bine folosite.

Începem cu câteva sugestii directe:

1.Particularitățile integrabile pot fi adeseori eliminate printr-o schimbare de variabilă. De exemplu , comportamentul singular K(t,s)s1/2 sau s-1/2 lângă s=0 poate fi eliminat prin transformarea z=s1/2. A se nota că presupunem că acest comportament singular este limitat la K , pe când cvadratura presupune de fapt produsul K(t,s)f(s) și tocmai acest produs trebuie să fie “fixat”. În mod ideal , în cel mai bun caz trebuie să deducem natura particulară a rezultatului înainte de a încerca o soluție numerică și să alegem acțiunea adecvată. De obicei un nucleu singular nu produce o soluție singulară f(t). ( Cel mai singular nucleu K(t,s)=(t-s) este pur și simplu operatorul de identitate , de explemplu ).

2.Dacă K(t,s) poate fi factorizat prin w(s)(t,s) , unde w(s) este singular și )(t,s) este omogenă, atunci cvadratura Gaussiană bazată pe w(s) ca o funcție de valori va merge bine. Chiar dacă factorizarea este aproximativă , convergența este adeseori îmbunătățită foarte mult. Tot ceea ce avem de făcut este să înlocuim gauleg în procedura fred 2 printr-o altă procedură de cvadratură. O secțiune anterioară a explicat cum să construim asemenea cvadraturi. De asemenea trebuie să-l înlocuim pe K cu .

Această metodă a un caz special al metodei produsului Nystrom[3,4] unde unii factori exteriori au termen singular p(t,s) depinzând și de t și de s din K și construiește valori potrivite pentru cvadratura sa Gaussiană. Calculele în cazul general sunt cam stânjenitoare , deoarece valorile depind de {ti} ales și de forma lui p(t,s).

Preferăm să implementăm rezultatul (produsul) Nystrom într-o grilă uniformă , cu o schemă de cvadratură care generalizează regula extinsă a lui Simpson la funcțiile de valoare oarecare (arbitrară). Vom discuta despre acestea în subsecțiunea de mai jos.

3.Formulele speciale de cvadratură sunt de asemenea folositoare când nucleul nu este strict singular , doar “aproape” singular. Un exemplu ar fi acela când nucleul e concentrat lângă t=s într-o proporție mult mai mică decât aceea în care soluția f(t) variază. În acest caz o formulă de cvadratură se bazează pe aproximarea într-un spațiu limitat al lui f(s) print-un polinom sau spline , în timp ce se calculează primele câteva “momente” ale nucleului K(t,s) la punctele de netezire ti .Într-o astfel de schemă dimensiunea mică a nucleului devine mai degrabă o calitate decât o responsabilitate. Cvadratura devine exactă când întinderea nucleului tinde spre zero.

4.O serie infinită de integrale e de asemenea o formă particulară. Reducerea seriei la o valoare finită mare ar trebui folosită doar ca o ultimă soluție. Dacă nucleul tinde rapid către zero , atunci o cvadratură Gauss-Laguerre [wexp(-s)] sau cvadratura Gauss-Hermite [wexp(-s2)] ar trebui să meargă bine.funcțiile lungi (stufoase) nu rezistă la transformarea

(3.1)

care duce pe 0<s< la 1>z>-1 astfel încât integrarea Gauss-legendre poate fi folosită. Aici >0 este o constantă pe care o putem ajusta pentru a îmbunătăți convergența.

5.O situație comună în practică este aceea că K(t,s) este singular de-a lungul diagonalei t=s. Aici metoda Nystrom dă greș deoarece nucleul este evaluat la (ti,si). Scăderea particularității este o soluție posibilă:

= (3.2)

unde r(t)= este calculat analitic sau numeric. Dacă primul termen din partea dreaptă este regulat putem folosi metoda Nystrom. În calculul ecuației (1.4) avem:

(3.3)

Uneori procesul de scădere trebuie repetat înainte ca nucleul să fie regular în întrgime.

Cvadratura pe o rețea uniformă cu pas variabil

În general este posibil să se găsească o cvadratură liniară de n puncte care să aproximeze integrala unei funcții f(x) , stabilește o funcție arbitrară cu valoarea w(x) , pe un interval oarecare (a,b) ca suma de valori să stabilească pe n chiar valori la distanțe egale ale funcției f(x) , să zicem la x=kh,(k+1)h, ….., (k+n-1)h. Schema generală pentru derivarea unei asemenea cvadraturi este să se noteze cele n ecuații liniare care trebuiesc rezolvate dacă cvadratura este exactă pentru cele n funcții f(x)=const,x,x2, ….., xn-1, și apoi să le rezolvăm pentru coeficienții respectivi. Acest lucru poate fi făcut analitic , o dată pentru totdeauna dacă momentele funcției de valoarea aceluiași interval de integrare

(3.4)

sunt cunoscute. Aici factorul h-n este ales pentru a face proporțional pe Wn cu h dacă (așa ca în cazul obișnuit) b-a este proporțional cu h.

Executând la bun sfârșit după metodă pentru un caz în patru puncte dă rezultatul:

(3.5)

În timp ce termenii din paranteză apar proporționali cu k2 există o neutralizare tipică și la O(k2) dar și la O(k).

Ecuația (3.5) poate fi specializată pe mai multe opțiuni ale lui (a,b)

Alegerea evidentă este a=kh, b=(k+3)h, în care caz obținem o cvadratură de 4 puncte care generalizează regula 3/8 a lui Simpson. De fapt putem recupera acest caz special dacă w(x)=1 , în care caz (3.4) devine

(3.6)

Cei patru termeni din ecuația cu paranteze drepte (3.5) devine fiecare independent de k și (3.5) se reduce de fapt la

(3.7)

Înapoi la cazul general al lui w(x) , alte alegeri pentru a și b sunt folositoare , de asemeni. De exemplu poate vrem să alegem ca (a,b) să fie ([k+1]h,[k+3]h) sau ([k+2]h,[k+3]h) , aceasta permițându-ne să terminăm o regulă extinsă al cărei număr de intervale nu este un multiplu de trei , fără să pierdem acuratețea: integrala va fi estimată folosindu-se cele patru valori f(kh),…, f([k+3]h). Chiar mai folositor ar fi să alegem ca (a,b) să fie ([k+1]h,[k+2]h), astfel folosind patru puncte pentru a integra un singur interval central. Aceste valori, atunci când sunt legate împreună într-o formulă extinsă dau scheme de cvadratură care au coeficienții (smooth) ,i.e. fără alternarea Simpson 2,4,2,4,2.

Toate aceste reguli sunt de același grad cu regula extinsă a lui Simpson , care este exactă pentru f(x) un polinom cubic. Regulile de grad inferior , dacă se dorește , se obțin în același mod similar. Formula de trei puncte este

(3.8)

Aici se ia cazul simplu w(x)=1 astfel încât

(3.9)

Atunci ecuația (3.8)devine regula lui Simpson

(3.10)

Pentru funcțiile de valoare variabilă w(x) , oricum ecuația (3.8) dă reguli cu un grad mai puțin decât regula lui Simpson , o dată ce nu beneficiază de simetria care este în cazul de valori constante.

Formula de două puncte este simplă

(3.11)

Aici există o procedură wwghts care folosește formulele de mai sus pentru a returna o cvadratură extinsă de N-puncte pentru intervalul (a,b) = (0, N-1îh). Datele de început wwghts sunt furnizate de utilizator prin procedura kermom apelată pentru a fi prima pentru momentele indefinite – integrale ale lui w(x) , anume

m=0,1,2,3 (3.12)

(limita inferioară este arbitrară și poate fi aleasă din conveniență). A se lua în considerare: când este apelat cu N<4 , wwghts returnează o regulă de rang inferior regulii lui Simpson ; trebuie să structurăm problema pentru a evita acest lucru.

int j,k;

double wold[5],wnew[5],w[5],hh,hi,c,fac,a,b;

hh=h;

hi=1.0/hh;

for (j=1;j!=n;j++) wghts[j]=0.0;

if (n >= 4){

b=0.0;

for (j=1;j!=n3;j++) {

c=j1;

a=b;

b=a+hh;

if (j == n3) b=(n1)*hh;

for (fac=1.0,k=1;k!=4;k++,fac*=hi)

w[k]=(wnew[k]wold[k])*fac;

wghts[j] += (

((c+1.0)*(c+2.0)*(c+3.0)*w[1]

(11.0+c*(12.0+c*3.0))*w[2]

+3.0*(c+2.0)*w[3]w[4])/6.0);

wghts[j+1] += (

(c*(c+2.0)*(c+3.0)*w[1]

+(6.0+c*(10.0+c*3.0))*w[2]

(3.0*c+5.0)*w[3]+w[4])*0.5);

wghts[j+2] += (

(c*(c+1.0)*(c+3.0)*w[1]

(3.0+c*(8.0+c*3.0))*w[2]

+(3.0*c+4.0)*w[3]w[4])*0.5);

wghts[j+3] += (

(c*(c+1.0)*(c+2.0)*w[1]

+(2.0+c*(6.0+c*3.0))*w[2]

3.0*(c+1.0)*w[3]+w[4])/6.0);

for (k=1;k!=4;k++) wold[k]=wnew[k];

}

}

else if (n == 3) {

(*kermom)(wnew,hh+hh,3);

w[1]=wnew[1]wold[1];

w[2]=hi*(wnew[2]wold[2]);

w[3]=hi*hi*(wnew[3]wold[3]);

wghts[1]=w[1]1.5*w[2]+0.5*w[3];

wghts[2]=2.0*w[2]w[3];

wghts[3]=0.5*(w[3]w[2]);

} else if (n == 2) {

wghts[1]=wnew[1]wold[1](wghts[2]=hi*(wnew[2]wold[2]));

}

}

Acum vom da un exemplu .

Exemplu : Un nucleu singular diagonal

Ca un caz particular , se consideră ecuația integrală

(3.13)

cu nucleul greoi (ales în mod arbitrar)

(3.14)

care are o singularitate logaritmică în stânga diagonalei , combinată cu o discontinuitate de rădăcină pătrată în dreapta.

Primul pas este să facem (în mod analitic în acest caz ) integralele cerute pentru fiecare indice pe partea singulară a nucleului , ecuația (3.12). o dată ce aceste integrale sunt date la o valoare fixă a lui x , putem folosi și pe x ca limită inferioară. Pentru orice valoare specificată a lui y , integrala nedefinită cerută este ori

(3.15)

sau

(3.16)

(unde o schimbare a variabilei a fost făcută în a doua egalitate din fiecare caz). Făcând aceste integrale în mod analitic (de fapt am folosit o grupare de integrare simbolică ) grupăm formulele rezultate în ordinea următoare. Observăm că w(j+1) returnează Fj(y;x).

#include <math.h>

extern double x;

void kermom(double w[], double y, int m)

Returnează în w[1..m] primii m indici de integrale nedefinite ale unei linii singulare a nucleului. (Pentru acest exemplu m este stabilit 4). Variabila z care va fi introdusă va reprezenta coloana , în timp ce variabila globală x va reprezenta linia. Îl putem alege pe x ca fiind limita inferioară a integrării. Astfel se va returna valoarea integralelor pentru fiecare m ori numai din partea stângă ,ori numai din partea dreaptă a diagonalei.

{

double d,df,clog,x2,x3,x4,y2;

if (y ?= x) {

d=yx;

df=2.0*sqrt(d)*d;

w[1]=df/3.0;

w[2]=df*(x/3.0+d/5.0);

w[3]=df*((x/3.0 + 0.4*d)*x + d*d/7.0);

w[4]=df*(((x/3.0 + 0.6*d)*x + 3.0*d*d/7.0)*x+d*d*d/9.0);

} else {

x3=(x2=x*x)*x;

x4=x2*x2;

y2=y*y;

d=xy;

w[1]=d*((clog=log(d))1.0);

w[2] = 0.25*(3.0*x+y2.0*clog*(x+y))*d;

w[3]=(11.0*x3+y*(6.0*x2+y*(3.0*x+2.0*y))

+6.0*clog*(x3y*y2))/18.0;

w[4]=(25.0*x4+y*(12.0*x3+y*(6.0*x2+y*

(4.0*x+3.0*y)))+12.0*clog*(x4(y2*y2)))/48.0;

}

}

Apoi scriem o procedură care construiește matricea de cvadratură.

#include !math.h?

#include ''nrutil.h''

#define PI 3.14159265

double x;

void quadmx(float **a, int n)

Construiește în a[1..n][1..n] matricea de cvadratură pentru un exemplu de ecuații Fredholm de tipul 2. Partea nesingulară a nucleului este calculată în această procedură în timp ce valorile de cvadratură care integrează partea singulară a nucleului sunt obținute în urma apelării lui wwghts. O procedură externă kermom care furnizează valorile integralei nedefinite la anumite momente ale părții singulare a nucleului trece la wwghts.

{

void kermom(double w[], double y, int m);

void wwghts(float wghts[], int n, float h,

void (*kermom)(double [], double ,int));

int j,k;

float h,*wt,xx,cx;

wt=vector(1,n);

h=PI/(n1);

for (j=1;j!=n;j++) {

x=xx=(j1)*h;

wwghts(wt,n,h,kermom);

cx=cos(xx);

for (k=1;k!=n;k++) a[j][k]=wt[k]*cx*cos((k1)*h);

++a[j][j];

}

free—vector(wt,1,n);

}

În sfârșit rezolvăm sistemul linear pentru orice parte dreaptă , aici sin x.

#include <stdio.h>

#include <math.h>

#include ''nrutil.h''

#define PI 3.14159265

#define N 40

int main(void) /* Program fredex */

Acest program dat ca exemplu arată cum să rezolvăm o ecuație Fredholm de tipul 2 folosind metoda produsului Nystrom și o cvadratură construită special pentru un anumit nucleu singular.

{

void lubksb(float **a, int n, int *indx, float b[]);

void ludcmp(float **a, int n, int *indx, float *d);

void quadmx(float **a, int n);

float **a,d,*g,x;

int *indx,j;

indx=ivector(1,N);

a=matrix(1,N,1,N);

g=vector(1,N);

quadmx(a,N);

ludcmp(a,N,indx,&d);

for (j=1;j!=N;j++) g[j]=sin((j1)*PI/(N1));

lubksb(a,N,indx,g);

for (j=1;j!=N;j++) {

x=(j1)*PI/(N1);

printf(''%6.2d %12.6f %12.6f„n'',j,x,g[j]);

}

free—vector(g,1,N);

free—matrix(a,1,N,1,N);

free—ivector(indx,1,N);

return 0;

}

Cu N=40 , acest program are o precizie de 10-5. Precizia crește când N4 (așa cum ar trebui pentru schema noastră de cvadratură de ordine Simpson) în ciuda mărimii nucleului singular. Figura (3.1) arată soluția obținută , de asemenea schițând soluția pentru valorile mai mici ale lui N, care sunt ele însele văzute ca fiind corecte. Se remarcă faptul că soluția este mai ușoară , chiar dacă nucleul este singular , ceea ce este un fenomen obișnuit.

CAPITOLUL 3

REZOLVAREA ECUAȚIILOR INTEGRALE

CU ALGORITMUL LUI RICHARDSON

3.1 ALGORITMUL LUI RICHARDSON

Fie Sk(n) polinom de interpolare care are valoarea xp în punctul tp , pentru p=n,n+1,…….,n+k, n,k N

Adică Sk(n) (tp) =xp , npn+k., gr Sk(n) k (nr.nodurilor –1)

Polinoamele Sk(n) ,n,k verifică următoarea relație de recurență

(Formula Aitken-Neville) (3.1)

Fie Rk(n) =Sk(n) (0)

Din (3.1) și t=0 rezultă

de unde rezultă

(3.2)

Această formulă (3.2) permite calculul numerelor ( prin recurență) Rk(n) , dacă se cunosc valorile inițiale R0(n) , nN. Dar, ,

polinom de grad 0 care ia valoarea xn în punctul t.

În particular R0(n) (0)=xn.

Așadar R0(n)=xn nN. (3.3)

Algoritmul descris de formulele descrise anterior , algoritm ce permite calculul șirurilor este cunoscut sub numele de Algoritmul lui Richardson.

Numărul Rk(n) se așează într-un tablou (tabloul Richardson) astfel

Înaintarea în acest tablou este de la stânga la dreapta și de sus în jos după așa numita regulă a triunghiului

Șirul (tn)n este un șir de numere distincte dat , iar șirul (xn)n este un șir a cărui convergență este accelerată de algoritmul lui Richardson.

Dacă în prima coloană a tabloului Richardson se cunosc termenii atunci ultimul termen din tablou care poate fi calculat este Rp(0).

Teorema1. Condiția necesară și suficientă pentru ca Rk(n) să fie x, nn0 este ca să existe numerele a1,a2,…….,akR a.î. .

Demonstrație. Fie – polinomul de interpolare cu proprietățile

Presupunem

Reciproc

Dacă

Dar și

teorema este demonstrată

Fie (tn)n un șir de numere reale convergent la zero, f o funcție indefinit derivabilă pe un interval I care conține șirul tn și xn=f(tn). Deoarece f este continuă șirul (xn)n este un șir convergent. Fie

Sk(n) fiind un polinom de interpolare al funcției f pe nodurile : tn, tn+1, …….,tn+k, avem

(3.4)

(formula restului în formule de interpolare)

Sk(n) polinom de interpolare al lui f pentru că

Pentru t=0 și (3.4) rezultă

(3.5)

Deoarece din (3.5)

De asemenea, întrucât șirul tn0 putem presupune ( eliminând un număr finit de termeni din acest șir) Atunci din (3.4) rezultă

Am demonstrat astfel următoarea teoremă

Teorema 2. Dacă xn=f(tn) , nN, unde f este o funcție indefinit derivabilă pe un interval I care conține șirul (tn)n și atunci și

Teorema 3. Fie (tn)nR+ , descrescător și convergent la zero.

.

Teorema 4. Fie (tn)nR+ , convergent la zero pentru care

Condiția necesară și suficientă pentru ca șirul (Rk+1(n))n să conveargă la mai repede decât șirul (Rk(n))n este ca :

Demonstrație

Din ipoteză

(3.6)

teorema este demonstrată

Algoritmul lui Richardson are și o generalizare

Fie crescătoare și continuă la dreapta în 0. Dacă (tn)nR+ , descrescător și converge la zero, atunci algoritmul generalizat al lui Richardson este descris de formulele

(3.7)

Forma clasică (3.2) și (3.3) rezultă di (3.7) pentru .

Rezultatele teoretice precedente se mențin înlocuind peste tot

STUDIUL COMPORTAMENTULUI SOLUȚIEI

APROXIMATIVE PENTRU O ECUAȚIE INTEGRALĂ

FREDLOLM DE SPEȚA A DOUA

Considerăm ecuația integrală următoare

(3.1)

Presupunem ca g să fie continuă și de două ori continuu diferențiabilă pe[0,1] , ca K(x,y) să fie continuă și să aibă derivate parțiale continue până la ordinul 2. Presupunem că (3.1) admite o soluție unică. Se știe atunci că această soluție admite o derivată secundară continuă.

Pentru a calcula integrala se folosește formula trapezelor cu un pas h=1/N

(3.2)

cu

În fine se vor da lui x valorile xi (i=0,1,2,…….,N) și se vor obține valori aproximative fi pentru f(xi) rezolvând sistemul liniar de N+1 ecuații cu N+1 necunoscute

(3.3)

Pentru a simplifica facem ipoteza

(3.4)

Scăzând (3.2) și (3.3) se obține pentru ei=fi-f(xi)

(3.5)

care este un sistem liniar de forma

Adaptând pentru vectori norma se obține pentru matricele pătratice corespunzătoare

este un vector de normă

hK este o matrice pătratică de normă (după (3.4)).

Dintr-o teoremă clasică se știe că I-hK admite un invers și că

Se deduce de aici că

Teorema 1. Dacă se presupune că, K și g au derivate secundare continue ( derivate parțiale pentru K) , că soluția lui (3.1) este unică , că, K verifică (3.4) , soluția aproximativă prin (3.3) există și este unică și eroarea verifică

unde K este o constantă independentă de I și h

Remarca

Ipoteza (3.4) a fost făcută pentru a simplifica demonstrația. Se poate în fond demonstra fără a utiliza (3.4) că pentru h suficient de mic

matricea I-hK este inversabilă și că norma inversei sale este mărginită în raport cu h (L.V.Kantorovich, Functional Analysis and Applied Mathematics, editată prin G.E. Forsythe).

Teorema 2. Cu aceleași ipoteze ca în teorema 1 înlocuim (3.3) prin

cu și C constantă. Avem atunci

( K* independent de i și de h) ; ( este suficient să se ia

).

Demonstrația este analoagă celei de la teorema 1.

Teorema 3. În plus față de ipotezele teoremei 1 se presupune acum existența și continuitatea derivatelor de ordinul 3 și 4 ale lui K și g ( derivate parțiale pentru K).

Eroarea verifică atunci

unde B este o funcție mărginită de i și de h și e(x) soluția ecuației integrale

(6.7)

Demonstrație: Soluția exactă f(x) verifică

cu

cu

Scăzând din (3.3) se obține punând :

ce corespunde rezolvării de la (3.7) prin metoda (3.6).

Se poate aplica teorema 2 și se obține

ceea ce demonstrează teorema 3.

Teorema 4. În plus față de ipotezele teoremei 3 presupunem existența și continuitatea derivatelor de ordinul 5 și 6 și g. Eroarea verifică atunci

unde A este o funcție mărginită de i și h, și (x) este soluția ecuației integrale

Demonstrația este comparabilă cu cea a teoremei 3.

Exemplu de extrapolare pentru o ecuație integrală

Considerăm ecuația integrală

Condiția (3.4) este verificată. În plus toate derivatele parțiale există până la un ordin oarecare. Căutăm soluția pentru un punct fix x=0.5 ; soluția exactă f(0.5) =1,006856727.

Soluția aproximativă w(h) obținută cu un pas h admite o desfășurare pară până la un ordin arbitrar

Se folosește procedeul Mr . Încercăm pași succesivi într-un raport 2 apoi într-un raport 1.5.

Prima experiență (=2)

A doua experiență (=1.5)

Se constată că coeficientul =1.5 conduce la aceeași precizie ca =2 dar cu un cost mai mic.

Pentru a atinge o precizie de 10-9 fără extrapolare (soluție normală) se va ajunge la rezolvarea unui sistem liniar de dimensiune excesivă.