Metoda Shooting

CUPRINS

1. Probleme la limită

1.1 Introducere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .4

1.2 Metoda shooting simplă . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .5

1.3 Metoda shooting simplă pentru probleme liniare la limită. . . . . . . 8

1.4 Teorema de Unicitate și Existență pentru soluția problemelor la limită . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .15

1.5 Dificultăți în executarea metodei shooting simple . . . . . . . . . . . . . .16

1.6 Metoda shooting multiplă . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

1.7 Indicații pentru implementarea metodei shooting multiple . . . . . 23

1.8 Un exemplu: Programul optim de control pentru ridicarea unei rachete spațiale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

1.9 Cazul al metodei shooting multiple ( Metoda lui Newton generalizata, Cvasiliniarizarea) . . . . . . . . . . . . . . . . . . . . . . . . . . . .30

2. Metode cu diferențe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3. Metode variaționale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

Bibliografie

CAPITOLUL 1

Probleme la limită

1.1 Introducere

Mult mai generale decât problemele cu valoare initială sunt problemele la limită. În cadrul acestora, se caută o soluție y(x) a sistemului de ecuații diferențiale ordinale cu n necunoscute,

(1.1.1a) y =f(x,y), y=, f(x,y)=,

satisfăcând o condiție de limită de forma :

(1.1.1b) Ay(a)+By(b)=c.

Aici, ab sunt numere date, A și B matrice pătratice de ordin n și c un vector din . În practică, condițiile de limită sunt de obicei separate :

(1.1.1) Ay(a)=c, By(b)=c,

adică, în (1.1.1b) liniile matricei [A,B,c] pot fi permutate astfel încât pentru matricea rearanjată [,,],

[,,]=.

Condițiile de limită sunt liniare în y(a), y(b).

Ocazional, în practică, se întâlnesc de altfel condițiile de limită neliniare de tipul

(1.1.1) r(y(a),y(b))=0,

care sunt formate cu ajutorul unui vector r de n funcții , i= ,de 2n variabile:

r(u,v)=.

Cu toate că problemele cu valoare ințială de obicei se rezolvă în mod unic, problemele la limită pot de asemenea să nu aibă soluție sau mai multe soluții.

Exemplu. Ecuația diferențială

(1.1.2a) w+w=0

pentru funcția reală , cu notația y(x)=w(x), y(x)=(x), poate fi scrisă sub forma (1.1.1a),

.

Are soluția generală

, c,c arbitrar alese.

Soluția specială w(x)=sinx este unica soluție care satisface condițiile de limită

(1.1.2b) w(0)=0, w()=1.

Toate funcțiile w(x)=csin x, cu c arbitrar, satisfac condițiile de limită

(1.2c) w(0)=0 , w()=0,

atâta timp cât nu există soluția w(x) din (1.2a) îndeplinind condițiile de limită

(1.2d) w(0)=0, w()=1.

(Se observă că toate condițile de limită (1.1.2b-d) au forma (1.1.1) cu A=B=[1,0].)

Multe probleme practice importante pot fi reduse la probleme la limită (1.1.1). Astfel este cazul, de exemplu, pentru problemele cu valori proprii pentru ecuații diferențiale, în care partea dreaptă a lui f a sistemului cu n ecuații diferențiale depinde de un parametru ,

(1.1.3a)

și trebuie satisfăcute n+1 condiții de limită de forma

(1.1.3b) r(y(a), y(b), A) = 0,

Problema (1.1.3) este determinată și de aceea, în general, nu are soluție pentru o alegere arbitrară a lui . Problema cu valori proprii din (1.1.3) constă în determinarea acelor numere, valorile proprii ale lui (1.3), pentru care (1.1.3) are o soluție. Prin introducerea unei funcții adiționale

și o ecuație diferențială adițională

,

(1.1.3) este echivalentă cu problema

, ,

care acum are forma (1.1.1), cu

, ,

mai mult decât atât, așa-numitele probleme la limită cu limită liberă sunt de asemenea reductibile cu (1.1). În aceste probleme, numai o singură abcisă a, este fixată, în timp ce b urmează sa fie determinat așa încât sistemul cu n ecuații diferențiale ordonate

(1.1.4a)

are o soluție y satisfacând cele n+1 condiții de limită

(1.1.4b) r(y(a),y(b))=0 ,

Aici, în locul lui x, se introduce o nouă variabilă independentă t și o constantă , care urmează să fie determinată, cu ajutorul

, ,

.

[Aici, în locul acestei alegerii a parametrilor, orice substituție de forma cu este de asemenea potrivită. Pentru , unde y(x) este o soluție a lui (1.4), se obține și (1.4) este astfel echivalent cu o problemă la limită de tipul (1.1) pentru funcțiile ,

(1.1.5) ,

, i=1,2,…,n+1

1.2 Metoda shooting simplă

Pentru început, dorim să explicăm metoda shooting simplă cu ajutorul unui exemplu.

Presupunem că ni se dă problema la limită

(1.2.1) ,

,

cu condiții de limită separată. Problema cu valoare inițială

(1.2.2) ,

în general are o soluție unic determinată w(x)=w(x;s) care bineînțeles depinde de alegerea inițială s pentru . Trebuie să se găsească a funcției F(s)=w(b;s)-. Pentru aceasta, trebuie determinată valoarea w(b)=w(b;s) a soluției w(x;s) a problemei cu valoare inițială (2.2) în punctul x=b. Calcularea lui F(s) este astfel echivalentă cu soluția unei probleme cu valoare inițială.

Deoarece w(b;s), și de aici F(s), sunt în generalle de limita pentru a de tipul ()pinde de parametrul 00000000000000000000000000000000000000 funcții diferențiabile continue a lui s, se poate folosi metoda lui Newton pentru a determina . Pornind cu o aproximație inițială atunci trebuie calculate iterativ valorile

astfel

(1.2.3)

Figura 14 Metoda shooting simplă

w(b;s), și astfel F(s), poate fi determinat rezolvând problema cu valoare inițială

(1.2.4) , , .

Valoarea derivatei lui F,

,

pentru poate fi obținut, de exemplu, alăturând o problemă cu valoare inițială suplimentară: cu ajutorul lui (719), se verifică ușor că funcția satisface

(1.2.5) , v(a)=0, .

Din cauza derivatelor parțiale , problema cu valoare inițială (1.2.5) este în general substanțial mai complexă decât (1.2.4). Din acest motiv, adesea se înlocuiește derivata în formula lui Newton (1.2.3), cu coeficientul de diferență ,

,

unde este ales destul de mic. este calculată, asemenea lui , rezolvând problema cu valoare initială. Pot apărea următoarele dificultăți :

Daca este ales prea mare, este o aproximație insuficientă pentru și

(1.2.3a) ,

converge către considerabil mai încet decât (1.2.3). Dacă este ales prea mic, atunci , și scăderea este supusă la anulare, astfel încât și cele mai mici erori din calcularea lui deteriorează puternic rezultatul

Figura 15 Graficul lui F(s)=w(1 ;s)-1.

Exemplu. Considerăm problema la limită

, w(0)=4, w(1)=1.

Urmărind (1.2.2), se găsește soluția problemei cu valoare inițială

, w(0 ;s)=4, .

Graficul lui F(s)=w(1;s)-1 este prezentat în Figura 15. Se observă că F(s) are două rădăcini . Iterația conform (2.3a) produce

Figura 16 prezintă graficele celor două soluții a problemei la limită. Soluțiile sunt calculate cu aproximativ 10 zecimale. Ambele soluții, întâmplător, pot fi exprimate în formă restrânsă de

,

Figura 16. Solutiile

unde reprezintă funcția eliptică iacobiană cu coeficient

.

Din teoria funcțiilor eliptice, pentru constantele se obțin valorile

Pentru soluția unei probleme la limită generală (1.1.1a), (1.1.1) care implică n funcții necunoscute

(1.2.6) r(y(a),y(b))=0, ,

unde f(x,y) și r(u,v) sunt vectori a n funcții, se procedează ca în exemplul de mai înainte. Se încearcă din nou determinarea unui vector de plecare pentru problema cu valoare inițială

(1.2.7) , y(a)=s

în așa fel că soluția y(x)=y(x;s) îndeplinește condițiile de limită ale lui (1.2.6),

r(y(a;s),y(b;s))r(s,y(b;s))=0.

Astfel trebuie gasită soluția a ecuației

(1.2.8) F(s)=0, F(s)r(s,y(b;s))=0

pe care o găsim cu ajutorul metodei lui Newton generală

(1.2.9) .

La fiecare pas al repetării, trebuie calculat , matricea iacobiană

DF()=,

și soluția a sistemului liniar de ecuații . Pentru calcularea lui trebuie determinat , adică, rezolvată problema cu valoare inițială (1.2.7) pentru . Pentru calcularea lui se observă

(1.2.10)

cu matricile

(1.2.11) , ,

.

În cazul funcțiilor neliniare r(u,v), oricum, nu se va calcula DF(s) cu ajutorul acestor formule complicate, dar în schimb se va aproxima cu ajutorul coeficienților de diferență. Astfel DF(s) va fi aproximat de matricea

unde

(1.2.12)

Ținând seama ca F(s)=r(s,y(b;s)), calcularea lui bineînțeles va solicita ca și să fie determinate prin soluția care corespunde problemelor cu valoare inițială.

Pentru condițiile de limită liniare (1.1.1b)

r(u,v)Au+Bv-c, ,

formulele se simplifică oarecum. Avem

F(s)As+By(b ;s)-c,

DF(s)A+BZ(b;s).

În acest caz, în locul formulei DF(s), este nevoie să determinăm matricea

.

Așa cum tocmai s-a descris, coloana j a lui Z(b;s) este înlocuită de coeficientul de diferență

.

Se obține aproximația

(1.2.13) , .

Deci, pentru a duce la bun sfârșit metoda lui Newton aproximativă

(1.2.14 )

trebuie făcut ceea ce urmează:

alegem un vector de pornire s.

Pentru i=0,1,2,…:

determinăm rezolvând problema cu valoare inițială (1.2.7) pentru , și calculăm F()=r()

alegem (suficient de mici) numerele și determinam rezolvând cele n probleme cu valori inițiale (1.2.7) pentru .

Calculăm cu ajutorul lui (1.2.12) [sau (1.2.13)] și de asemenea soluția a sistemului liniar de ecuații , și punem .

La fiecare pas al metodei trebuie astfel rezolvate cele n+1 probleme cu valoare inițială și al n-lea sistem de ecuații liniare.

Ținând seama de convergența locală simplă a metodei lui Newton (aproximativă) (1.2.14), metoda va diverge în general la o soluție a lui F(s)=0, dacă nu vectorul de pornire este deja suficient de aproape de o soluție a lui F(s)=0.

1.3. Metoda shooting simplă pentru problemele la limită liniare

Substituind cu din (1.2.9), se pierde în general convergența locală de gradul II a metodei lui Newton. Metoda de substituție (1.2.14), de obicei, converge numai liniar (local) rata convergenței fiind mai mare pentru aproximații mai bune ale lui cu .

În cazul special al problemelor cu valoare de limită liniară, avem acum pentru orice s (și pentru alegerea arbitrară a lui ), astfel încât (1.2.9) și (1.2.14) devin identice. Printr-o problemă liniară de limită se înțelege o problemă în care f(x,y) este o funcție afină în y și condițiile de limită (1.1.1b) sunt liniare,

(1.3.1)

Ay(a)+By(b)=c,

cu o matrice de T(x), o funcție , și matricile constante A și B. Presupunem în continuare ca T(x) șiproximația

(1.2.13) , .

Deci, pentru a duce la bun sfârșit metoda lui Newton aproximativă

(1.2.14 )

trebuie făcut ceea ce urmează:

alegem un vector de pornire s.

Pentru i=0,1,2,…:

determinăm rezolvând problema cu valoare inițială (1.2.7) pentru , și calculăm F()=r()

alegem (suficient de mici) numerele și determinam rezolvând cele n probleme cu valori inițiale (1.2.7) pentru .

Calculăm cu ajutorul lui (1.2.12) [sau (1.2.13)] și de asemenea soluția a sistemului liniar de ecuații , și punem .

La fiecare pas al metodei trebuie astfel rezolvate cele n+1 probleme cu valoare inițială și al n-lea sistem de ecuații liniare.

Ținând seama de convergența locală simplă a metodei lui Newton (aproximativă) (1.2.14), metoda va diverge în general la o soluție a lui F(s)=0, dacă nu vectorul de pornire este deja suficient de aproape de o soluție a lui F(s)=0.

1.3. Metoda shooting simplă pentru problemele la limită liniare

Substituind cu din (1.2.9), se pierde în general convergența locală de gradul II a metodei lui Newton. Metoda de substituție (1.2.14), de obicei, converge numai liniar (local) rata convergenței fiind mai mare pentru aproximații mai bune ale lui cu .

În cazul special al problemelor cu valoare de limită liniară, avem acum pentru orice s (și pentru alegerea arbitrară a lui ), astfel încât (1.2.9) și (1.2.14) devin identice. Printr-o problemă liniară de limită se înțelege o problemă în care f(x,y) este o funcție afină în y și condițiile de limită (1.1.1b) sunt liniare,

(1.3.1)

Ay(a)+By(b)=c,

cu o matrice de T(x), o funcție , și matricile constante A și B. Presupunem în continuare ca T(x) și g(x) sunt funcții continue pe [a,b]. Din y(x;s) din nou deducem soluția problemei cu valoare ințială

(1.3.2) y(a;s)=s

Pentru y(x;s) se poate da o formulă explicită,

(1.3.3) y(x ;s)=Y(x)s+y(x ;0),

unde matricea de Y(x) este soluția problemei cu valoare inițială

, Y(a)=I

Din partea dreaptă a lui (3.3) deducem

u(a ;s)=Y(a)s+y(a ;0)=Is+0=s,

adică, (x;s) este o soluție a lui (1.3.2). Deoarece, sub ipotezele pentru T(x) și g(x) făcute mai sus, problema cu valoare inițială are o soluție unică, rezultă că u(x;s)=y(x;s). Folosind (1.3.3), se obține pentru funcția F(s) din (1.2.8)

(1.3.4) F(s)=As+By(b;s)-c=[A+BY(b)]s+By(b;0)-c.

Astfel, F(s) este de asemenea o funcție afină a lui s. Prin urmare [conform (1.2.13)],

.

Soluția a lui F(s)=0 [presupunând existența lui ] este dată de

sau, oarecum mai general, de

(1.3.5)

unde este arbitrar. Cu alte cuvinte, soluția a lui F(s)=0) va fi aflată prin metoda (2.14), printr-un singur pas al iterației, inițiat cu un vector de pornire .

1.4 Teorema de Unicitate și Existență pentru soluția problemelor la limită

erminam matricea

corespunde 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

Sub condiții foarte restrictive se poate arăta singura rezolvare a anumitor probleme la limită. În final, consideram următoarele probleme la limită cu condiții de limită neliniare:

(1.4.1) ,

Problema dată în (4.1) este rezolvată cu precizie dacă funcția F(s) din (1.2.8) are o rădăcină :

(1.4.2)

Aceasta din urmă este cu siguranță adevarată dacă se poate găsi o matrice nesingulară Q de astfel încât

(1.4.3)

este o aplicație contractivă din ; rădăcina a lui F(s) este un punct fixat a lui, ()=.

(1.4.4) Teoremă. Pentru problemele la limită (1.4.1), următoarele ipoteze sunt satisfacute:

(1) f si sunt continue pe

(2) Există un kC[a,b] cu pentru orice (x,y)s.

(3) Matricea

P(u,v)=r(u ,)+ r(u,)

admite pentru orice u,v o reprezentare de forma

P(u,v)=P(I+M(u,v))

cu o matrice constantă nesingulară și o matrice M=M(u,v), și există constantele și m cu

pentru orice u,v .

(4) Există un număr cu astfel încât

Atunci problema la limită (1.4.1) are exact o soluție y(x).

Demonstrație. Arătam că dacă facem o alegere potrivită a lui Q, și anume Q=, funcția din (1.4.3) satisface

(1.4.5) pentru orice ,K=

Din acesta rezultă imediat că:

pentru orice ,

este o aplicație contractivă care, are exact un punct fix , care este o rădăcină a lui F(s), în raport cu nesingularitatea lui Q.

Pentru avem

(1.4.6)

=

=

=,

unde matricea

este soluția problemei cu valoare inițială

, Z(a;s)=I,

În funcție de ipotezei (2), astfel rezultă că urmează pentru Z

estimarea

și mai departe, din (4.6) și ipotezele (3) și (4),

=

Teorema este acum demonstrată.

Observație. Condițiile teoremei sunt numai suficiente si sunt de asemenea foarte restrictive. Ipoteza (3), de exemplu, deja in cazul n=2, nu este satisfacuta pentru astfel de conditiile de limita simple ca .

1.5 Dificultăți în executarea metodei shooting simple

De la metoda de rezolvare a problemei la limită

r(y(a),y(b))=0

se va cere ca pentru fiecare din intervalul de definiție a soluției y se va produce o valoare aproximativă pentru . În metoda shooting discutată până acum, numai valoarea inițială va fi determinată. Problema este astfel rezolvată, deoarece valoarea a soluției la fiecare punct poate fi aproximativ determinată considerând problema cu valoare inițială

(1.5.1) ,

Aceasta, oricum, este adevarată numai teoretic. În practică, greșelile cresc considerabil dacă soluția y(x)=y(x;) a lui (1.5.1) depinde de . Un exemplu va demonstra aceasta.

Exemplu1. Sistemul liniar de ecuații diferențiale

(1.5.2)

are soluția generală

(1.5.3) , arbitrare,

Fie y(x;s) soluția lui (1.5.2) satisfăcând condiția inițială

.

Se verifică imediat că

(1.5.4)

Dorim să determinăm soluția y(x) a lui (7342) care satisface condițiile liniare de limită separate

(1.5.5) y(0)+ y(10)-=,

sau mai simplu,

, .

Soluția exactă se deduce din (1.5.3)

(1.5.6)

Valoarea inițială a soluției exacte este

.

În calcularea lui , de exemplu, în 10 puncte aritmetice, în locul lui se obtine cea mai bună valoare aproximativă a formulei

Cu.fie, . La valoarea inițială

Oricum, se găsește conform (1.5.4) o soluție exactă cu

.

În raport cu (1.5.4) avem în acest exemplu pentru x>0 suficient de mare,

,

influența datelor inexacte inițiale crește în mod exponențial cu x. Pentru soluția y(x;s) a problemei cu valoare ințială avem

,

Această estimare, oricum, ne arată de asemenea că influența datelor inițiale inexacte s=y(a) poate fi arbitrar redusă printr-o reducere a intervalului . Situația dificilă ar fi astfel eliminată daca s-ar cunoaște valorile , k=1,2,…,m, a soluției y(x) la mai multe puncte suficient de aproape cu

Alegând astfel încât să nu fie prea mare, se poate determina y() pentru fiecare , adică la o exactitate suficientă rezolvând problema cu valoare inițială

pe un interval mic cu metoda din sectiunea 7.2. Astfel apare problema cum poate fi calculata valoarea pentru un particular. Pentru aceasta, ca pentru valoarea inițială , se poate aplica din nou metoda shooting simplă. Întra-adevăr, deducând din soluția problemei cu valoare inițială

(1.5.7)

este o rădăcină a funcției

.

Această rădăcină din nou poate fi determinată iterativ (, ), de exemplu cu ajutorul metodei lui Newton modificată. Singura diferență față de metoda shooting simplă constă în faptul că fiecare pas din metoda lui Newton pentru determinarea lui , trebuie calculate două soluții, și anume și , cu ajutorul soluțiilor problemei cu valoare inițială de tipul (1.5.7). Pentru valoarea sugerează că valoarea de pornire pentru metoda iterativă Newton, pe când teoretic (cu condiția ca este egal cu și este calculat fără a rotunji erorile). Presupunând însă că este suficient de mic, efectul erorilor în pe , este ținut în limite, .

Aceasta, oricum, nu îndepărtează toate dificultățile. În calcularea lui , k=1,2,…,m, problema următoare, de fapt, frecvent apărută, care de mai multe ori restricționează semnificația practică a metodei shooting simplă: funcția f din ecuația diferențială , pentru a fi siguri, uneori are derivate parțiale continue în raport cu y pentru , dar este nelimitat pe . În acest caz, soluția a problemei cu valoare inițială , este nevoie să fie definită numai într-o vecinătate oarecare a lui , a cărei lungime depinde de s. Astfel, si posibil să existe numai pentru valorile s într-o mulțime mică M.

Exemplul 2. Considerăm problema la limită [conform Troesch(1960), (1976)]

(1.5.8) ,

(1.5.9) y(0)=0 , y(1)=1

( un parametru fixat ).

Ținând cont că problema poate fi tratată cu metoda simplă de tragere, trebuie prima data ‘estimată‘ panta inițială . În integrarea numerică a problemei cu valoare inițială (1.5.8) cu, se arată că soluția y(x ;s) depinde extrem de mult de s: pentru s=0.1,0.2,…calculul se termină chiar înainte ca limita de dreapta (x=1) să fie atinsă, trebuie ca exponentul să crească, y(x ;s) are un singur punct (care depinde de s) la un Efectul pantei inițiale pe amplasarea particularității (locația singularității) poate fi estimată în exemplul următor :

contine prima integrală

(1.5.10) .

Condițiile definesc constanta de integrare,

.

Integrarea lui (1.5.10) conduce la

.

Atunci punctul singular este dat de

.

Pentru o evaluare aproximativă a integralei descompunem intervalul de integrare,

arbitral ,

și estimăm integralele parțiale separat. Găsim

și,

astfel se obține estimarea

pentru fiecare cantitatea este o limită superioară pentru ; deci în particular,

pentru orice .

Comportarea asimtotică a lui pentru poate fi găsită cu usurință: Pentru , avem în prima aproximație

astfel încât, asimtotic pentru ,

(1.5.11)

[Se poate chiar arăta (vezi mai jos) că, asimtotic pentru ,

(1.5.12)

se poate realiza. Dacă, în procesul shooting, se ajunge la limita de dreapta x=1, valoarea lui poate fi astfel aleasă cel mult atât de mare că

i.e. ;

pentru se obține un domeniu mic

(1.5.13)

Pentru ‘expert’ adăugăm următoarele: Problema cu valoare inițială

(1.5.14) ,

are soluția exactă

;

aici sn și cn sunt funcții iacobiene eliptice cu modulele k, care aici depinde de panta inițială s. Daca K(k) indică perioada împătrită a lui cn, atunci cn are o rădăcină (1.5.15)

și de aceea y(x ;s) are o particularitate logaritmică acolo. Dar K(k) are o extindere

sau, rescris ținând cont de s,

din care rezultă (1.5.12).

Pentru soluția problemei la limită actuală,adică,

y(1 ;s)=1,

se găsește pentru valoarea

s = 4.575.04614*10

[conform (5.13)].Această valoare este obținută dintr-o relație dintre funcțiile iacobiene combinate cu o metodă iterativă. Pentru completare, adăugăm încă o relație în plus. Soluția problemei la limită (1.4.8), (1.4.9) are o particularitate logaritmică

(1.5.16)

pentru rezultă că

1.0326… ;

particularitatea soluției exacte se întinde în vecinatatea imediată al punctului de limită din partea dreaptă. Pentru o soluție numerică directă a problemei fără cunoașterea teoriei funcțiilor eliptice, vezi secțiunea 1.7.

1.6 Metoda shooting multiplă

Metoda shooting multiplă a fost descrisă de mai multe ori în literatură, de exemplu în Keller (1968), Osborne(1969), și Bulirsch (1971). Un program FORTRAN poate fi găsit în Oberle și Grimm (1989).

Într-o metodă multiplă shooting, valorile

k =1…,m,

a soluției exacte y(x) a problemei la limită

(1.6.1)

în mai multe puncte

sunt calculate simultan prin iterații. Fie y(x ;x,s) soluția problemei cu valoare inițială

Problema constă acum în determinarea vectorilor s, k=1,2,…,m, în așa fel că funcția pentru , , ,

pusă împreună cu , este continuă și astfel o soluție a ecuației diferențiale , și în plus satisface condițiile de limită r(y(a),y(b))=0 (vezi figura 17). Aceasta produce următoarele condiții nm :

(1.6.2) , ,

pentru componentele necunoscute nm a lui :

Cu totul , (1.6.2) reprezintă un sistem de ecuații de forma:

(1.6.3) ==0

în necunoscutele :

Figura 17. Metoda shooting

Se poate rezolva iterativ cu ajutorul metodei lui Newton,

(1.6.4)

La fiecare pas al metodei trebuie calculat F(s) și DF(s) pentru . Pentru calcularea lui F(s) trebuie determinate valorile pentru k=1,2,…,m-1 rezolvând problemele cu valoare inițială

,

și să calculăm F(s). Matricea iacobiană

,

ținând seama de structura specială a lui are forma

(1.6.5)

unde matricile A, B, , de , k=1,…,m-1, pe rând sunt matrice iacobiene,

(1.6.6) k=1,2,…,m-1,

Este eficient în practică a înlocui coeficienții diferențiali din matricele A,B, cu coeficienții de diferență care pot fi calculați rezolvând cele (m-1)n probleme cu valoare inițială (n probleme cu valoare inițială pentru fiecare matrice G,…,). Calcularea lui din conform lui (6.4) poate fi realizată după cum urmează : Cu notațiile

(1.6.7) ,

(1.6.4) este echivalentă cu următorul sistem linear de ecuații:

(1.6.8)

Începând cu prima ecuație, se pot exprima toate în funcție de . Astfel se găsește

(1.6.9)

și din acestea în final, cu ajutorul ultimei ecuații

(1.6.10)

unde

Acesta este un sistem de ecuații liniare pentru un vector necunoscut , care poate fi rezolvat cu ajutorul eliminării lui Gauss. Odată determinat , se obține ,…, succesiv din (6.8) și din (1.6.7). Mai departe remarcăm că sub ipotezele Teoremei (4.4) se poate arata că matricea din (6.10) este nesingulară. Va fi din nou o matrice nesingulară Q de nmnm astfel încât funcția

contractează, și astfel are exact un punc fix cu F=0.

Este arătat, în plus că F(s), și cu atât mai mult DF(s) este definită pentru toți vectorii

astfel încât iterația (1.6.4) a metodei shooting multiple poate fi executată pentru . Aici , k=1,2,…,m-1, este mulțimea tuturor vectorilor , pentru care soluția există (cel puțin ) pe intervalul de dimensiune mică . Această mulțime include toate mulțimile al tuturor pentru care există pe toate intervalele [a,b]. Metoda lui Newton pentru calcularea lui , cu ajutorul metodei shooting simple poate fi executată numai pentru . Aceasta arată că cererile metodei shooting multiple pentru valoarea vectorilor de pornire în metoda lui Newton sunt considerabili mai simple decât cele de la metoda shooting simplă.

1.7. Indicatii pentru implementarea a metodei shooting multiple

Metoda shooting multiplă este mai degrabă costisitoare. Iterația în metoda lui Newton modificată , de exemplu, are forma

(1.7.1) , ,

și la fiecare pas trebuie calculată cel puțin aproximația a matricei iacobiene DF(s) formând coeficienți de diferență corespunzători. În continuare considerăm problema alegerii punctelor intermediare , k=1,2,…,m, din [a,b], cu condiția ca o soluție aproximativă (traiectoria de pornire) este cunoscută pentru problema la limită. Punem =a . având deja ales (<b), integrăm problema cu valoare inițială

,

cu ajutorul metodelor din secțiunea 7.2, și terminăm integrarea în primul punct pentru care soluția devine “prea mare ”comparativ cu , adică 2; apoi avem .

Considerăm problema la limită

(1.7.2) , .

Pentru traiectoria de pornire tragem o linie dreaptă unind cele două puncte de limită [].

Subdiviziunea găsită de calculator este aratată în figura 18. Începând cu această subdiviziune și valorile metoda shooting multiplă produce soluția cu 9 cifre în 7 iterații.

Problema sugerează o traiectorie de pornire ‘‘mai avantajoasă’’ : problema liniarizată

are soluția . Această funcție ar fi trebuit să aibă o alfel de traiectorie de pornire ‘‘mai bună’’.

Figura 18. Subdiviziunea în metoda shooting multiplă

În unele probleme practice partea dreaptă a funcției f(x,y) a ecuației diferențiale este numai pe porțiuni continuă pe [a,b] ca o funcție de parametru x, sau altfel continuă pe porțiuni diferențiabile. În aceste cazuri trebuie să includem punctele de discontinuitate printre subdiviziunea punctelor ; altfel sunt probleme la convergență (vezi Exemplul 2).

Astfel rămâne problema cum să găsești o primă aproximație pentru problema la limită. În multe cazuri se cunoaște comportarea calitativă a soluției astfel încât se poate găsi ușor cel puțin o aproximație brută. De obicei și aceasta este de ajuns pentru convergență că metoda lui Newton modificată nu impune în special cereri înalte pe calitatea vectorului de pornire.

În cazurile complicate așa-numitele metode de continuare (metode de omotopii) sunt de obicei eficace. Aici se merge treptat pe dibuite, aproape toate problemele conțin parametri siguri ,

(1.7.3)  :

și problema constă în determinarea soluției y(x), care bineînțeles depinde de asemenea de , pentru o valoare sigură . Este în general adevărat că pentru o altă valoare sigură problema la limită este simplă, asfel încât cel puțin o aproximație bună este cunoscută pentru soluția a lui . Metoda de continuare constă acum în următoarele : începând cu , prima dată se determină soluția a lui . Apoi se alege o secvență finită de numere suficient de apropiate cu , , și pentru avem soluția a lui ca aproximație de pornire pentru determinarea soluției a lui . Apoi în în final avem soluția dorită a lui . Ca metoda să meargă este important să nu alegem pași prea mari , și ca parametri naturali , care sunt intrinseci în problemă, să fie aleși astfel: introducerea parametrilor artificiali de tipul , unde f(x,y) este dat de partea dreaptă a ecuației diferențiale și g(x,y) o funcție arbitrar aleasă simplă care nu are nici o legatură cu problema, nu dă rezultate în cazurile critice.

O variantă mult mai precisă a metodei de continuare, care s-a dovedit a fi utilă în contextul metodelor shooting multiple, este descrisă în Deuflhard (1979).

Exemplul2. (O problemă singulară la limită). Temperatura de distribuție din interiorul unui cilindru de raza 1 este descrisă de soluția a problemei neliniare la limită

(1.7.4)

aici, este un parametru natural cu

, .

Pentru calcul numeric al soluției pentru se poate folosi subdiviziunea

și, începând cu soluția cunoscută explicită construim lanțul de omotetie pentru y(x ;0.8). Problema, oricum expune o mai mare dificultate: partea dreaptă a ecuației diferențiale are o particularitate la x=0. Atâta timp cât este adevarat că și soluția este definită, într-adevăr analitic, pe tot intervalul , se ridică totuși dificultăți de convergență considerabile. Motivul este descris în următoarele: o analiză înapoi arată că rezultatul al calcului numeric poate fi interpretat ca soluție exactă a problemei la limită (1.7.4) cu datele de limită ușor modificate

Această soluție este apropiată de soluția , dar spre deosebire de are numai o derivată continuă în x=0, deoarece

Printr-un artificiu dificultățile pot fi evitate. Soluțiile învecinate cu sunt alese cu ajutorul unei dezvoltări a seriei de puteri în punctul x=0:

(1.7.5) ,

Coeficienții , i=2,3,4,…, pot fi toți exprimați în funcție de , o constantă necunoscută; prin substituția în ecuația diferențială (1.7.4) se găsește

(1.7.6)

din care, pentru ,

Diferențierea dă

cu

Mai departe,

și

Se poate arăta că , și în general, .

Problema singulară la limită poate fi tratată acum după cum urmează. Pentru se folosește în vecinătatea lui x=0 reprezentarea seriei de puteri (1.7.5), și la o distanță suficientă de o particularitate x=0 ecuația diferențială (1.7.4),

(1.7.7)

Acum partea dreaptă încă conține parametrul necunoscut . Așa cum este arătat în Secțiunea 1), aceasta poate fi interpretată ca o problemă la limită extinsă. Se notează:

Figura 19.Solutia pentru .

și se obține un sistem de ecuații diferențiale pentru ,

(1.7.8)

cu condițiile de limită

(1.7.9)

În metoda shooting multiplă se alege punctul neuniform ca unul din subdiviziunea de puncte

Nu este nici un motiv să ne ingrijorăm în legatură cu planeitatea soluției în punctul neuniform; este asigurată prin construcție. Cu ajutorul omotetiei [vezi (1.7.3)], începând cu , soluția este ușor de obținut din din numai 3 iterații cu o exactitate de aproape 6 cifre (timpul total de calcul pe un calculator CDC 3600 pentru toate 8 traiectorii : 40secunde). Rezultatul este arătat în figura 19.

1.8.Un exemplu: Programul optim de control pentru ridicarea unei rachete spațiale

Următoarea problemă extinsă își are originea în călătoria în spațiu. Acest exemplu concludent ilustrează cum problemele neliniare la limită sunt mânuite și rezolvate în practică, deoarece experiența arată că soluția actuală ale acestor probleme reale cauzează cele mai mari dificultăți pentru un începător. Soluția numerică a fost dusă la bun sfârșit cu metoda shooting multiplă; un program poate fi găsit în Oberle și Grimm (1989).

Coordonatele traiectoriei ale unei rachete de tipul Apollo satisfac, în timpul unui zbor al rachetei prin admosfera pământului, următoarele ecuații diferențiale :

(1.8.1)

Notațiile sunt:-viteza; -unghiul direcției de zbor; h-altitudinea de deasupra suprafeței pământului ; R-raza pământului ; -altitudinea normalizată; – distanța de la suprafața pământului, -densitatea admosferică; -modulul de rezistență aerodinamic; -coeficientul forței portante aerodinamică; u-parametrul de control, care poate fi ales arbitrar ca o funcție de timp; g-accelerația gravitațională; S/m -(suprafața frontală)/(masa rachetei).

Valorile numerice sunt  : ;  ;;

;(vezi figura 20).

Ecuațiile diferențiale au fost simplificate oarecum presupunând (1) un pământ sferic nemișcat, (2) o traiectorie de zbor într-un mare plan de cerc, (3) ca racheta și astronauții pot fi supuși la o deaccelerare nelimitată. Partea dreaptă a ecuațiilor diferențiale conține totuși termeni care fundamental sunt fizici. Cel mai mare efect este produs de termenii multiplicați de și respectiv ; aceștia sunt forțe admosferice care, în ciuda densității mici a aerului , sunt semnificanți ținând seama de viteza mare a rachetei; ei pot fi influentați prin parametrul u (unghiul critic). Termenii multiplicați de g sunt forțe gravitaționale ale pământului care acționează asupra rachetei; restul de termeni rezultă din alegerea sistemului de coordonate.

Figura 20.Coordonatele traiectoriei

În timpul zborului prin atmosfera pământului, racheta este încinsă considerabil. Punctul mort total convectiv încălzit pe unitate de suprafață este dat de integrala

(1.8.2)

Rata integrării este de la t=0, timpul când racheta atinge nivelul de 400,00ft admosferice, până la un timp instant T. Racheta urmează să fie manevrată într-o poziție inițială favorabilă pentru amerizarea finală în Pacific. Prin parametrul liber disponibil u manevrarea urmează se să execute astfel încât, căldura J devine minimală și ca următoarele condiții de limită sunt satisfacute: Datele în momentul intrării:

(1.8.3)

Datele la sfârșitul manevrării :

(1.8.4)

Timpul terminal T este liber. este ignorat în procesul de optimizare.

Calculul variațional [cf,eg,Hestenes(1966)] ne învață următoarele: se formează, cu parametrii (multiplicatorii Lagrange), expresia :

(1.8.5)

unde satisfac cele trei ecuații diferențiale

(1.8.6)

Controlul optim u este dat de

(1.8.7)

Se observă că (1.8.6), ținând cont de (1.8.7), este neliniar în .

Deoarece timpul terminal T nu este expus la nici o condiție, mai departe condițiile de limită trebuie satisfăcute :

(1.8.8) .

Astfel problema este redusă la o problemă la limită pentru cele 6 ecuații diferențiale (1.8.1), (1.8.6) cu cele 7 condiții de limită (1.8.3), (1.8.4), (1.8.8). Astfel avem de a-face cu o problemă la limită liberă [cf (1.1.4a,b)]. O soluție în formă apropiată este imposibilă; trebuie folosite metode numerice.

Lipsa de experiență nu ar trebui să fie indusă în eroare de forma inocentă a părții drepte a ecuației diferențiale (1.8.1): în timpul integrării numerice se observă rapid că depind foarte mult de datele inițiale.

Construirea unei traiectorii de pornire. Din motive aerodinamice graficul parametrului de control u va avea forma reprezentată în Figura 21 (informație de la inginerii spațiului). Această funcție poate fi aproximată, de exemplu, de

(1.8.9)

unde

și , pentru moment sunt constante necunoscute.

Figura 21. Parametrul de control

Pentru determinarea lui se rezolvă următoarele probleme la limită auxiliare: ecuațiile diferențiale (1.8.1) cu forma u din (1.8.9) și în plus

(1.8.10)

Cu condițiile de limită (1.8.3), (1.8.4) și T=230.

Valoarea de limită T=230 secunde este o aproximare pentru durata manevrei de reîntoarcere. Cu aproximațiile relativ insuficiente

Problema la limită auxiliară poate fi rezolvată chiar și cu metoda shooting simplă cu condiția să se integreze în direcția înapoi (punctul inițial , punctul terminal ). Rezultatul după 11 iteratii este

(8.11)

Soluția problemei la limită actuală. Cu funcția de control ’neoptimă’ din (8.9), (8.11) se obțin aproximații la integrând (8.1). Această traiectorie de pornire ‘incompletă’ poate fi făcută completă dupa cum urmează: deoarece cos u >0, rezultă din (1.8.7) că ; alegem

.

În continuare, din

rezultă o aproximație pentru . O aproximație pentru poate fi obținută din relația , unde H este dat de (1.8.5), deoarece H=constant este o primă integrală a ecuațiilor mișcării.

Figura 22.Traiectoriile h,

Figura 23 Traiectoriile ale variabilelor alaturate

Figura 24.Aproximatiile succesive pentru controlul u

Cu această traiectorie aproximativă pentru (T=230) și tehnicile din Secțiunea 1.7, se pot determina acum punctele de subdiviziune pentru metoda shooting multiplă (m=6 ajunge).

Figura 22 și Figura 23 arată rezultatul calculării (dupa 14 iterații); timpul total al manevrării de întoarcere optimă pare să fie P=224.9 sec; precizia rezultatelor este de 7 cifre.

Figura 24 arată comportamentul controlului în timpul desfășurării iterațiilor. Poate fi văzut cum săriturile inițial mari din punctele de subdiviziune ale metodei shooting multiple sunt ‘netezite’.

Cazul a metodei shooting multiple

(Metoda lui Newton generalizat[, Cvasiliniarizarea)

Dacă subdiviziunea intervalului [a,b] este făcută foarte bine (), metoda multiplă shooting converge către metoda lui Newton generală pentru problemele la limită [conform, de exemplu, Collatz(1966), de altfel cunoscută ca cvasiliniarizarea].

În această metodă, o aproximație pentru soluția exactă y(x) a problemei neliniare la limită (6.1), este îmbunătățită rezolvând problemele liniare la limită. Folosind deyvoltarea Taylor a funcției f(x,y(x)) și r(y(a),y(b)) pentru ,se găsește în prima aproximație

unde Se așteaptă, prin urmare ca soluția a problemei liniare la limită

(1.9.1)

să fie o soluție mai bună a lui (6.1) decât . Pentru utilizarea întâziată introducem funcția de corecție care prin definiție este soluția problemei cu valoare inițială

(1.9.2)

asfel încât

(1.9.3)

Înlocuind din 1.9.1 cu , se obține o aproximație în plus pentru y. În ciuda simplei derivări, această metodă de iterație are serioase dezavantaje care aduc câteva îndoieli în legătură cu valoarea ei practică:

Funcțiile vectorului ,…trebuie să fie conținute peste întreg intervalul lor [a,b].

Funcția matrice trebuie calculată într-o formă explicit analitică și trebuie de asemenea să fie conținută peste întreg intervalul ei [a,b]. Realizarea ambelor cereri este aproape imposibilă în problemele care apar în practică simultan

Dorim să arătăm că metoda multiplă shooting cu converge la metoda 1.9.1 astfel: fie o funcție pe [a,b] și fie mai departe (1.9.1) unic rezolvabilă cu soluția Apoi următoarele sunt adevărate:

dacă în metoda multiplă shooting (1.6.3)-(1.6.10) se alege orice subdiviziune a lui și vectorul de pornire

cu

atunci soluția lui (1.6.8) va produce un vector

(1.9.4) cu .

Pentru demonstrație presupunem pentru simplicitate Dorim să arătăm că pentru fiecare h există o funcție diferențiabilă asfel încât și

Prima dată arătăm că producțiile

ce apar în (1.6.9), (1.6.10) au limite simple cu (hk=const, hj= const). Reprezentarea (6.9) pentru poate fi scrisă în forma

(1.9.5)

conform lui (1.6.3), este dat de

astfel

și cu ajutorul teoremei de valoare medie,

(1.9.6)

Fie acum Z(x) soluția problemei cu valoare inițială

(1.9.7)

și matricea soluția lui

.

Asfel încât

(1.9.8)

pentru matricile se arată ușor prin inducție dupa k că

(1.9.9)

Într-adevăr, deoarece aceasta este adevarată pentru k=1. Dacă ipoteza este adevarată pentru k-1, atunci funcția

satisface ecuația diferențială și condițiile inițiale . Din unicitatea soluțiilor problemelor cu valoare inițială rezultă că și de aici (1.9.9).

Avem matricile

(1.9.10)

cu notația

prin scăderea lui (1.9.8) și (1.9.10) se obține aproximarea

(1.9.11)

Dacă pentru orice este Lipschitziană uniform continuă în y, se poate arăta usor pentru că

precum și mărginirea uniformă a lui pentru din (1.9.11).

– în particular, pentru

(1.9.12)

cu o constantă independentă de k și h.

Din identitatea

și aproximațiile următoare

,

independent de k și h, astfel obținem, tinând cont de și (1.9.9),

(1.9.13)

cu o constantă D independentă de k și h.

Din (1.9.5), (1.9.6), (1.9.9), (1.9.12), (1.9.13) se obține

(1.9.14)

Aici este funcția

(1.9.15)

Precizăm că depinde de asemenea de h, deoarece depinde de h. Evident, este diferențiabilă în x , și din (1.9.7) avem

astfel încât

(x)=.

Substrăgând această ecuație din (1.9.3), obținem pentru diferența

,

ecuația

(1.9.16)

și mai departe cu ajutorul lui (1.9.7)

(1.9.17) pentru

cu o alegere potrivită a lui K.

Ținând cont de (1.9.14) și (1.9.17) este suficient să arătăm că în loc să completăm demonstrația lui (1.9.4). Cu ajutorul lui (1.9.15) se arată în același fel ca (1.9.14) că satisface ecuația de forma

.

Din rezultă că

Scăderea (1.9.2), ținând seama de (1.9.17), produce

și astfel , deoarece (1.9.1) din ipoteză este unic rezolvabilă și de aici A+BZ(b) este nesingulară. Ținând cont de (1.9.14), (1.9.17), se demonstrează (1.9.4).

Capitolul 2

Metode cu diferențe

Ideea fundamentală de la baza tuturor metodelor cu diferențe este să înlocuim coeficienții diferențiali dintr-o ecuație diferențială cu coeficienții de diferență potriviți și să rezolvăm ecuațiile diferite care se obțin astfel.

Vom ilustra aceasta cu următoarea problemă la limită simplă de ordin ordinul II pentru o funcție

(2.1)

Sub ipotezele că q,g (q și g sunt funcții continue pe [a,b] și q(x) pentru ) poate fi arătat că (2.1) are soluția unică y(x).

Împărțim intervalul [a,b] în n+1 intervale egale,

,

si, cu notațiile înlocuim coeficientul diferențial pentru i= cu al doilea coeficient de diferență 

Vom estima eroarea Presupunem că funcția y este de 4 ori continuă diferențiabilă pe [a,b], . Apoi , din deyvoltarea lui Taylor a functiei , de coeficient, se găsește

și astfel

Deoarece este încă continuă, rezultă că

(2.2) ……….

Datorită lui (2.1), valorile satisfac ecuațiile

(2.3) i = 1,2,…,n

Cu notațiile vectorii

, ,

și matricea diagonală simetrică de

(2.4) A=,

ecuațiile (2.3) sunt echivalente cu

(2.5)

Metodele cu diferențe constau acum în scăderea termenului de eroare din (2.5) și luând soluția a sistemului de ecuații liniare rezultat,

(2.6) Au=k

ca o aproximație la .

Prima dată dorim să arătăm câteva proprietăți ale matricii A din (2.4). Vom scrie AB pentru două matrici de dacă pentru Avem următoarele

(2.7) Teoremă. Dacă pentru I=1,..,n, atunci A din (2.4) este pozitiv definită, și cu o matrice pozitiv definită de nn

(2.8)

Demonstrație. Începem arătând că A este pozitiv definită. Pentru aceasta, considerăm matricea de nn . Avem pentru valorile proprii aproximarea , și de aici . Dacă sunt valori proprii ale lui , rezultă că det; oricum, det, se verifică usor cu ajutorul formulei de recurență că det (care se obține imediat extinzând determinantul dupa primul rând). Valoarea nu este o valoare proprie a lui , deoarece ar fi singulară. Dar

D=diag(1,-1,1,…,),

astfel încât , și cu , de asemenea este nesingulară. Astfel obținem estimarea pentru valorile proprii ale lui . Aceasta arată în particular că este pozitiv definită. Din

pentru , si ,

rezultă imediat că pentru ; de aici definirea pozitivă a lui A. Rămâne să demonstrăm inegalitatea . În final,consideram matricile D, , D, J, cu

,

, .

Deoarece evident avem inegalitățile

(2.9)

Ținând seama de și de estimarea pentru valorile proprii ale lui , avem pentru valorile proprii a lui , pentru raza de acțiune a lui . Din

și deoarece , avem de asemenea convergența din

Din (2.9) , , avem

așa cum trebuia arătat.

Din teorema de mai sus rezultă în particular că sistemul de ecuații (2.6) are o soluție (dacă pentru ), care poate fi ușor găsit, cu ajutorul metodei lui Cholesky. Deoarece A este o matrice diagonală, numărul de operații pentru soluția lui (2.6) este proporțional cu n. Vrem acum să derivam o aproximație pentru erorile

ale aproximațiilor obținute din (2.6) pentru soluțiile exacte

(2.10) Teoremă. Fie problema la limită (2.1) care are o soluție și fie pentru . De asemenea fie pentru și fie soluția a lui (2.6). Apoi pentru i=1,2,..,n,

Demonstrație. Ținând cond de (2.5) și (2.6) avem pentru ecuația

Folosind ecuațiile

pentru

obținem din teorema (2.7) și reprezentarea (2.2) a lui ,

(2.11)

unde e=[1,1…1]. Vectorul poate fi obținut imediat cu observația următoare: problema la limită specială

, y(a)=y(b)=0,

de tipul (2.1) are soluția exactă . Pentru această problemă la limită, oricum avem din (2.2) și soluția discretă u din (2.6) coincide cu soluția exactă din (2.5). În plus pentru această problemă la limită specială matricea A din (2.4) este tocmai matricea din (2.8), și pe deasupra k=e. Astfel avem , . Împreună cu (2.11) se realizează ipoteza teoremei.

pag. 222-227

Capitolul 3

Metode variaționale

Metodele variaționale (metodele Rayleigh-Ritz-Galerkin) sunt bazate pe faptul că soluțiile unor importante tipuri de probleme la limită posedă cu siguranță proprietățiile minimale . Dorim să explicăm aceste metode pentru următoarea problemă la limită simplă pentru o funcție

(3.1)

, .

Se observă că problema (3.1) este într-un fel mai generală decât (2.1).

Sub ipotezele

(3.2) , ,

, ,

,

cu cea mai mică valoare proprie a problemelor cu valori proprii

z(a)=z(b)=0,

este cunoscut că (3.1) întotdeauna are exact o soluție. Pentru următoarele presupunem (3.2) și facem ipoteza simplificată g(x,u(x))=g(x).

Daca u(x) este soluția lui (3.1), atunci y(x) = u(x)-l(x) cu

, ,,

este soluția problemei la limită de forma 

(3.3)

y(a)=0, y(b)=0,

cu valori de limită care dispar. Fără pierderea generalității putem asfel considera, în locul lui (3.1), problemele de forma (3.3). Cu ajutorul operatorului diferențial

(3.4)

asociat cu (3.3), vrem sa formulăm problema (3.3) cumva diferit. Operatorul L trasează mulțimea

tuturor funcțiilor reale care sunt de două ori continue diferențiabile pe [a,b] și satisfac condițiile de limită în mulțimea C[a,b] a funcțiilor continue pe [a,b]. Problema la limită (3.3) este astfel echivalentă cu găsirea unei soluții a lui

(3.5) L(y)=f,

evident este un spațiu vectorial real și L un operator liniar pe : pentru de asemenea aparțin lui , și avem pentru toate numerele reale Pe mulțimea al funcțiilor pătrat integrabile pe [a,b] introducem acum o formă biliniară și o normă cu ajutorul definiției

(3.6) .

Operatorul diferențial L din (3.4) are câteva proprietăți care sunt importante pentru întelegerea metodelor variaționale. Avem următoarele:

(3.7) Teoremă. L este un operator simetric pe D, avem

pentru orice

Demonstrație. Prin integrarea în părti se găsește

,

deoarece u(a)=u(b)=0 pentru . Din motive de simetrie rezultă de asemenea că

(3.8)

și de aici afirmația.

Partea dreaptă a lui (3.8) nu este definită numai pentru .

Într-adevăr fie mulțimea tuturor funcțiilor u care sunt absolut continue pe intervalul [a,b] cu u(a)=u(b)=0, pentru care din [a,b] (există aproape peste tot) este de două ori integrabil.

În particular, toate funcțiile diferețiabile continue pe porțiuni care satisfac condițiile de limită aparțin lui D. D este un spațiu vectorial real cu Partea dreaptă a lui (3.8) definește pe D o formulă biliniară simetrică

(3.9)

care pentru coincide cu . Ca mai înainte, se arată pentru y, u, prin integrarea în părți, că

(3.10) .

Relativ la produsul scalar introdus din (3.6), L este un operator definit pozitiv în următorul sens :

(3.11).Teoremă. Sub ipotezele (3.2) avem

[u,u]=(u,L(u))>0 pentru orice .

Avem chiar estimarea

(3.12) ,pentru orice

cu norma și constantele

, .

Demonstrație. Ținând cont că > 0, este suficient să arătam (3.12). Pentru , deoarece u(a)=0,

pentru .

Inegalitatea lui Schwarz produce aproximația

și astfel

(3.13) .

Conform ipotezei (3.2), avem , pentru  ; din (3.9) și (3.13) rezultă astfel că :

.

În final din (3.13) avem de asemenea

,

așa cum trebuia arătat.

În particular, din (3.11) deducem imediat unicitatea soluției y din (3.3) sau (3.5). Dacă , atunci =0 și de aici , care produce imediat .

Definim acum pentru o funcțională patratică

(3.14) F(u)=[u,u]-2(u,f) ;

unde f este partea dreaptă a lui (3.3) sau (3.5). F asociază cu fiecare funcție un număr real F(u). Fundamental pentru metodele variaționale este observația că funcția F atinge cea mai mică valoare exactă pentru soluția y a lui (3.5) :

(3.15)Teoremă. Fie soluția lui (3.5). Atunci

F(u)>F(y)

Pentru orice u

Demonstrație. Avem L(y) = f și de aici din (3.10) și din definiția lui F, pentru

F(u)=[u,u]-2(u,f)=[u,u]-2(u,L(y))

=[u,u]-2[u,y]+[y,y]-[y,y]

=[u-y,u-y]-[y,y]

>-[y,y]=F(y),

deoarece, din teorema (3.11), [u-y,u-y] >0 pentru .

Ca un rezultat lateral, observăm identitatea

(3.16) [u-y,u-y]=F(u)+[y,y] pentru orice .

Aproximarea minimă a lui F poate fi obținută sistematic după cum urmează : se alege un subspațiu dimensional finit S a lui D, . Dacă din S = m, atunci referindu-ne la o bază a lui S, fiecare admite o reprezentare de forma

(3.17) , .

Apoi se determină minimul a lui F pe S,

(3.18)

Și se ia ca fiind o aproximație pentru soluția exactă y a lui (3.5), care conform lui (3.15) minimizează F pe tot spațiul D. Pentru calcularea aproximației considerăm următoarea reprezentare a lui F(u), , obținută prin (3.17),

.

cu ajutorul vectorilor și a matricei A de mm,

(3.19) , , ,

se obține pentru funcția pătratică

(3.20) .

matricea A este pozitiv definită deoarece A este simetrică din (3.9), și pentru toți vectorii avem de asemenea și astfel, din teorema (3.11),

.

sistemul de ecuații liniare

(3.21)

Prin urmare, are o soluție unică , care poate fi calculată cu ajutorul metodei Cholesky. Ținând cont de identitatea

și pentru , rezultă imediat că pentru , și rezultă în consecință că funcția

aparținând lui furnizează minimul (3.18) a lui F(u) pe S. Cu y soluția lui (3.5), rezultă imediat din (3.16), în funcție de F(u)=minF(u), că

(3.22)

Vrem să folosim această relație pentru a estima eroarea . Avem următoarele :

(3.23) Teoremă. Fie y soluția exactă a lui (3.3),(3.5). Fie un subspatiu finit dimensional a lui D, și fie . Atunci estimarea

merge pentru orice . Aici C =unde sunt constantele teoremei (3.11). Demonstrație.Pentru arbitral, avem

.

Din aceasta rezultă ipoteza.

Fiecare majorant superior pentru , sau mai ușor pentru , imediat apare o estimare pentru . Dorim să stabilim o asfel de limită pentru un caz special imporatant. Alegem pentru subspatiul S a lui D multimea

a tuturor funcțiilor spline cubice care aparțin unei subdiviziuni fixe a intervalului [a,b],

,

și eliminăm la sfârșit punctele a și b. Evident, . Notam cu lungimea celui mai mare subinterval al subdiviziunii ,

,

și mai departe mulțimea

(pentru subdiviziuni echidistante, K=1).

Funcția spline cu

, i=0,1,…,n,

pentru

unde y este soluția exactă a lui (3.3), (3.5), clar aparțin lui . Se obține estimarea

cu condiția ca . Împreuna cu teorema (3.23), aceasta dă următorul rezultat:

(3.24) Teoremă. Fie soluția exactă y a lui (3.3), aparține lui , și fie ipotezele (3.2), satisfăcute. Fie , și funcții spline pentru care

Atunci există constantele K și C care pot fi determinate independent de y astfel încât

Estimarea poate fi îmbunătățită :

eroarea de limită merge astfel la 0 ca putere terță a fineții lui ; în această privință, metoda este superioară față de metoda cu diferențe a secțiuni precedente [vezi teorema 3.10].

Pentru implementarea metodei variaționale, în cazul , prima dată trebuie selectată o bază pentru . Se vede ușor că m=dim [funcția spline este determinată în mod unic de cele n + 1 condiții

, I=1,2,…,n-1,

Pentru arbitrare. Se poate găsi usor o bază a funcțiilor spline din astfel încât

(3.25) pentru si .

Această bază are avantajul că matricea corespunzatoare A din (3.19), este o matrice de bandă de forma

(3.26)

deoarece în funcție de (3.9), (3.25),

oricând . Odată ce componentele matricei și ale vectorului din (3.19) pentru această bază au fost găsite prin integrare, se rezolvă sistemul liniar de ecuații (3.21) pentru și astfel se obține aproximația pentru soluția exactă y. Bineînteles, în locul lui se pot alege de asemenea alte spații De exemplu se poate lua pentru S mulțimea

al tuturor polinoamelor de grad cel mult n unde a și b dispar. Matricea A în cazul acesta, nu ar mai fi o matrice de bandă.

Similar cu funcțiile spline, se poate de asemenea alege pentru S mulțimea

,

unde este din nou o subdiviziune si este alcătuit din toate funcțiile care în fiecare subinterval ,I=0,1,…,n-1, coincide cu un polinom de grad mai 2m-1 (“spațiul funcției Hermite”). Aici din nou printr-o alegere potrivită a bazei {}din , se poate asigura că matricea A = ([]) este o matrice de bandă. Printr-o alegere potrivită a parametrului m se poate chiar obține metode de ordin mai mare decat 3 în felul acesta [conform teoremei (3.24)] pentru , analog cu teorema (3.24), se obține

,

cu condiția ca ; avem o metodă de ordin cel puțin 2m-1.

Cea mai mare parte a muncii din metodele de tipul acesta constă în calcularea coeficienților sistemului liniar de ecuații (3.21) cu ajutorul integrării. Pentru problemele la limită din ecuațiile diferențiale ordinale aceste metode de obicei sunt prea costisitoare și prea inexacte și nu concurează cu metoda multiplă shooting (vezi secțiunea următoare pentru rezultate comparative). Valorile ei devin evidente numai în problemele la limită pentru ecuații diferențiale parțiale.

Spațiile dimensional finite a funcțiilor ce satisfac condițiile de limită ale problemelor la limită (3.5) de asemenea joacă un rol în metodele de colocație. Ideea acestor metode este foarte naturală: se încearcă aproximarea soluției y(x) a lui (3.5) cu o funcție u(x) din S reprezentată de în funcție de o bază a lui S. În final, se selectează m puncte de colacatie diferite ,j=1,2,…,m și se determină astfel încât

(3.27a) , j=1,2,…m,

adică, ecuația diferențială L(u)=f este satisfacută exact în punctele de colocație. Desigur, aceasta este echivalentă cu rezolvarea următoarelor ecuații liniare pentru coeficienții ,k=1,2…,m:

(3.27b) , j=1,2…,m

Sunt multe posibilități de a implementa astfel de metode, și anume, prin alegerile diferite ale lui S, ale bazelor din S, și a punctelor de colocație . Cu o alegere potrivită, se pot obține metode foarte eficiente.

De exemplu, este avantajos pentru funcțiile bazelor ,j=1,2,…,m să avem un suport compact, cum ar fi B-splinele [vezi (3.25)]. Apoi matricea a sistemului liniar (3.27) este o matrice de bandă. Apoi, cu o alegere potrivită a punctelor de colocație , adesea se poate rezolva (3.27) cu ajutorul transformărilor rapide ale lui Fourier.

Pentru ilustrare,considerăm problema la limită simplă

de forma (3.3). Aici, toate funcțiile satisfac condițiile de limită. Dacă alegem punctele de colocație ,j=1,2,…,m, atunci din (3.27) numerele vor rezolva problema cu interpolare geometrică

, j=1,2,…,m.

Prin urmare, poate fi calculat folosind transformările rapide Fourier.

Bibliografie

Differential algebraic equation index transformations. SIAM J. Sci. Statist.

Comput. 9, 39-47 (1988).

Petzold, L. R.: ODE methods for the solution of differential/algebraic systems.

SIAM J. Numer. Anal. 21, 716-728 (1984).

Gottlieb, D., Orszag, S. A.: Numerical Analysis of Spectral Methods: Theory and

Applications. Philadelphia: SIAM (1977).

Gragg, W.: Repeated extrapolation to the limit in the numerical solution of ordinary

differential equations. Thesis, UCLA (1963).

: On extrapolation algorithms for ordinary initial value problems. J. SIAM

Numer. Anal. Ser. B 2,384-403 (1965).

Griepentrog, E., Marz, R.: Differential-Algebraic Equations and Their Numerical

Treatment. Leipzig: Teubner (1986).

Grigorieff, R. .D.: Numerik gewohnlicher Differentialgleichungen 1. 2. Stuttgart: Teub

ner (1972, 1977).

Hairer, E., Lubich, Ch.: Asymptotic expansions of the global error of fixed step-size

methods. Numer. Math. 45, 345-360 (1984).

, , Roche, M.: The numerical solution of differential-algebraic systems

by Runge-Kutta methods. In: Lecture Notes in Mathematics 1409. Berlin, Heidel

berg, New York: Springer (1989).

,Nersett, S. P., Wanner, G.: Solving Ordinary Differential Equations 1. NonstijJ

Problems. Berlin, Heidelberg, New York: Springer (1987).

– -, Wanner, G.: Solving Ordinary Differential Equations II. Stiff and Differential

Algebraic Problems. Berlin, Heidelberg, New York: Springer (1991).

Henrici, P.: Discrete Variable Methods in Ordinary Differential Equations. New York:

John Wiley (1962).

Hestenes, M. R.: Calculus of Variations and Optimal Control Theory. New York: John

Wiley (1966).

Heun, K.: Neue Methode zur Approximativen Integration der Differentialgleichungen

einer unabhangigen Variablen. Z. Math. Phys. 45, 23-38 (1900).

Horneber, E. H.: Simulation elektrischer Schaltungen auf dem Rechner. Fachberichte

Simulation, Bd. 5. Berlin, Heidelberg, New York: Springer (1985).

Hull, E. E., Enright, W. H., Fellen, B. M., Sedgwick, A. E.: Comparing numerical

methods for ordinary differential equations. SIAM J. Numer. Anal. 9, 603-637

(1972). [Errata, ibid. 11, 681 (1974).]

Isaacson, E., Keller, H. 8.: Analysis of Numerical Methods. New York: John Wiley

(1966).

Kaps, P., Rentrop, P.: Generalized Runge-Kutta methods of order four with stepsize

control for stiff ordinary differential equations. Numer. Math. 33, 55-68 (1979).

Keller, H. B.: Numerical Methods for Two-Point Boundary-Value Problems. London:

Blaisdell (1968).

Krogh, F. T.: Changing step size in the integration of differential equations using

modified divided differences. In: Proceedings of the Conference on the Numerical

Solution of Ordinary Differential Equations, 22-71. Lecture Notes in Mathematics

362. Berlin, Heidelberg, New York: Springer (1974).

Kutta, W.: Beitrag zur naherungsweisen Integration totaler Differentialgleichungen.

Z. Math. Phys. 46,435-453 (1901).

Lambert, J; D;: Computational Methods in Ordinary Differential Equations. London,

New York, Sydney, Toronto: John Wiley (1973).

Na, T. Y., Tang, S. c.: A method for the solution of conduction heat transfer with

non-linear heat generation. Z. Angew. Math. Mech. 49,45-52 (1969).

Oberle, H. 1., Grimm, W.: BNDsco-A program for the numerical solution of optimal

control problems. Internal Report No. 515-89/22, Institute for Flight Systems

Dynamics, DLR, Oberpfaffenhofen, Germany (1989)..

Bank, R. E., Bulirsch, R., Merten, K.: Mathematical modelling and simulation of electrical circuits and semiconductor devices. ISN M 93, Basel: Birkhauser ( 1990).

Boyd, J. P..: Chebyshev and. Fourier Spectral Methods. Berlin, Heidelberg, New York:

Springer (1989).

Bulirsch, R.: Die Mehrzielmethode zur numerischen Losung von nichtlinearen Rand

wertproblemen und Aufgaben der optimalen Steuerung. Report of the Carl

Cranz-Gesellschaft (1971).

, Stoer, J.: Numerical treatment of ordinary differential equations by extrapola

tion methods. Numer. Math. 8, 1-13 (1966).

Butcher, J. c.: On Runge-Kutta processes of high order. J. Austral. Math. Soc. 4,

179-194 (1964).

Byrne, D. G., Hindmarsh, A. c.: Stiff o.d.e.-solvers: A review of current and coming

attractions. J. Compo Phys. 70, 1-62 (1987).

Ciarlet, P. G., Lions, J. L. (Eds.): Handbook of Numerical Analysis. Vol. 11. Finite

Element Methods (Part 1). Amsterdam: North Holland (1991).

-, Schultz, M. H., Varga, R. S.: Numerical methods of high order accuracy for

nonlinear boundary value problems. Numer. Math. 9,394-430 (1967).

'–, Wagschal, c.: Multipoint Taylor formulas and applications to the finite

element method. Numer. Math. 17,84-100(1971).

Canuto, C, Hussaini, M. Y., Quarteroni, A., Zang, T. A.: Spectral Methods in Fluid

Dynamics. Berlin, Heidelberg, New York: Springer (1987).

Coddington, E. A., Levinson, N.: Theory of Ordinary Differential Equations. New

York: McGra~-Hi\\ {1955). Conatz, L.: The Numerical Treatment of di.fferential Equations. Berlin: Springer (1960). –: Functional Analysis and Numerical Mathematics. New York: Academic Press

(1966).

O:.~lq~~.;t 0.: [Jon.vergen.ce a.n.d ..ta.b;{~ty {n t(..e numerical integration of ordinary

differential equations. Math. Scand. 4, 33-53 (1956).

: Stability and error bounds in the numerical integration of ordinary differen

tial equations. Trans. Roy. Inst. Tech. (Stockholm), No. 130 (1959).

-: A special stability problem forlinearmultistep methods. 81T3, 27-43 (1963).