Curs 10.

Metode numerice de rezolvare a ecuaţiilor diferenţiale şi a sistemelor de ecuaţii diferenţiale

cu aplicaţii în ingineria electrică

METODE NUMERICE

Modelul matematic cel mai des întâlnit al fenomenelor care stau la baza majorităţii aplicaţiilor electrotehnice este ecuaţia diferenţială. Rezolvarea exactă a ecuaţiilor diferenţiale ordinare este posibilă pentru o clasă foarte restrânsă de ecuaţii!!!

O ecuaţie diferenţială este o ecuaţie care conţine pe lângă variabilele independente şi funcţiile necunoscute şi derivatele acestor funcţii (sau diferenţialele lor) până la ordinul n inclusiv (numărul n reprezintă ordinul ecuaţiei diferenţiale).

O ecuaţie diferenţială se numeşte ordinară dacă conţine o singură variabilă independentă şi are forma generală:

( ) 0y,...,''y,'y,y,xf )n( =

Comportarea dinamică a sistemelor fizice conduce la modele matematice formate din ecuaţii diferenţiale ordinare sau sisteme de ecuaţii diferenţiale care nu pot fi rezolvate pe cale analitică (funcţii complicate ca formă sau funcţii cunoscute doar pe baza unor valori în puncte date tabelar şi obţinute pe cale experimentală). Din acest motiv se recurge la rezolvarea numerică a acestora.

Metodele numerice de aproximare a soluţiilor conduc la tabele de valori ale funcţiei necunoscute. Valorile tabelate se calculează utilizând o valoare deja calculată cu un pas înainte (metode unipas) sau câteva valori calculate deja (metode multipas – metode de tip Adams).

Rezolvarea unei ecuaţii diferenţiale asociate unui circuit electric de ordin I excitat cu un impuls.

Condensator de capacitate C care se încarcă de la o sursă de tensiune continuă E, printr-un rezistor de rezistenţă R.

Descărcare unui condensator de capacitate C, încărcat iniţial la tensiunea E, pe un rezistor de rezistenţă R.

Analiza stabilităţii la mari perturbaţii a unui generator sincron legat la un sistem de putere infinită

Se consideră un sistem electroenergetic (SEE) de configuraţie simplă: generator sincron (GS), legat la un sistem de putere infinită printr-o reţea radială (transformator bloc şi linie dublu circuit de înaltă tensiune).

Enunţ:

Cunoscând parametrii elementelor de sistem şi datele referitoare la un anumit regim de funcţionare, se cere să se elaboreze un program de calcul pentru analiza stabilităţii la mari perturbaţii a GS prin integrarea numerică a ecuaţiilor diferenţiale care descriu funcţionarea în regim tranzitoriu a SEE.

Se vor considera numai perturbaţii simple, iar pentru integrarea sistemului de ecuaţii diferenţiale se vor utiliza metode monopas şi multipas, cu algoritm explicit sau implicit.

Modelul matematic

Analiza stabilităţii la mari perturbaţii se face prin integrarea ecuaţiei diferenţiale de mişcare a ansamblului rotoarelor generatorului şi turbinei:

dd t M

P P Hm e

2

21δ ω= ⋅ − − ⋅( )

rezultând curba de variaţie în timp a unghiului intern δ al generatorului (curba de oscilaţie) şi cea a vitezei unghiulare ω(reprezentând, de fapt, abaterea vitezei unghiulare faţă de turaţia sincronă ωs = 314 rad/s)

Analiza formei acestor curbe oferă informaţii în privinţa stabilităţii sau a instabilităţii generatorului la perturbaţia considerată.

M - constanta mecanică a ansamblului turbină-generator, Pm - puterea mecanică a GS, Pe - puterea electrică a GS, H - constanta de amortizare (înglobând efectele tuturor surselor de amortizare a oscilaţiilor).

dd tdd t M

P P Hm e

δ ω

ω ω

=

= ⋅ − − ⋅

⎧

⎨⎪⎪

⎩⎪⎪

1 ( )

Ecuaţia diferenţială se poate scrie sub forma unui sistem de două ecuaţii diferenţiale de ordinul întâi:

ttt

0

0 0

0 0

0=

=

=

⎧

⎨⎪⎪

⎩⎪⎪

ω ωδ δ

( )( )

unde variabila independentă este timpul t, iar puterea electrică are expresia:

P P Pe e e= + ⋅ −1 2 sin ( )δ α

unde Pe1, Pe2 şi α sunt constante, având valorile dependente de parametrii elementelor de sistemşi de regimul de funcţionare a sistemului!!!

Condiţiile iniţiale cunoscute sunt definite de relaţia

În concluzie, analiza stabilităţii tranzitorii se face prin rezolvarea numerică a sistemului de ecuaţii diferenţiale cu condiţiile iniţiale pentru t∈[ 0;tfinal ] şi interpretarea rezultatelor obţinute.

0.30

0.34

0.38

0 0.5 1 1.5 2 2.5t [s]

δ [rad]

-0.20

-0.10

0.00

0.10

0.20

0.30

0 0.5 1 1.5 2 2.5

t [s]

ω [rad/s]

Curba de variaţie în timp a vitezei unghiulare ω

Curba de variaţie în timp a unghiului intern δ

Metode numerice directe de rezolvare a ecuaţiilor diferenţiale

Soluţia m în punctul xi ( [a, b], a = x0< x1<…<xn= b) în cadrul acestei metode este calculata pe baza informaţiei conţinută numai în nodul precedent xi-1 şi eventual din intervalul [xi-1, xi].

Metoda dezvoltării în serie Taylor

Fie funcţia f:IxR cu valori in R, unde I este un interval real, f fiind o funcţie continua dată iar y0 fiind valoarea initiala.

Ne propunem să găsim funcţia y care satisface problema cu valori (condiţii) iniţiale (care se numeşte problema Cauchy):

⎩⎨⎧

==

00 y)x(y)y,x(f'y

bx...xxa n21 <<<<<

care constă în evaluarea funcţiei y(x) în nodurile aparţinând intervalului de definiţie

2=n ))(2

(01 yx fffhfhyy ⋅+⋅+⋅+=

⇒ n2 y,...,y

))(2

(1 yxii fffhfhyy ⋅+⋅+⋅+=+

Pentru

Prin recurenţă

Aproximaţia este cu atât mai bună cu cât numărul de termeni luaţi în considerare în dezvoltarea Taylor este mai mare. Metoda este directă întrucât pentru calculul lui yi+1sunt necesare informaţii numai despre punctul anterior (xi, yi).

xff x ∂

∂=

yff y ∂

∂=

Demonstratia 1. Pe tablă!!

Metoda lui Euler ( metoda clasică )

Este cea mai simplă metodă de integrare numerică a ecuaţiilor diferenţiale ordinare. Se obţine din metoda Taylor pentru n=1, adică se reţin numai primii doi termeni din dezvoltare rezultând forma explicită a metodei lui Euler:

),(1 iiii yxfhyy ⋅+=+ ,...1,0=i

nxb

h 0−=

Interpretare geometrică: se alege un pas de integrare h astfel încât intervalul de definiţie[x0, b] să fie împărţit în paşi egali:

)y,x(f'y = 00 y)x(y =

)x(yy =

Avem aceeaşi problemă de rezolvare a ecuaţiilor diferenţiale cu condiţii iniţiale:

,

şi avem şi curba soluţiei

( )ξ''2

2

yhE ⋅≤ 1+<< ii xx ξ

Prin metoda lui Euler soluţia în nodul xi+1 se aproximează cu ordonata punctului de intersecţie a tangentei la curbă în punctul (xi, yi) cu dreapta x=xi+1.

Ecuaţia tangentei:)(')( xyxxyy ii ⋅−+=

),()(' ii yxfxy =

rezultă formula de recurenţă a algoritmului Euler:),(1 iiii yxfhyy ⋅+=+

Metoda lui Euler se numeşte şi metoda liniilor poligonale pentru că curba y=y(x) se înlocuieşte prin linia poligonală M0, M1,… conform figurii 2 (dreapta ce trece prin M0 cu coeficientul unghiularf(x0,y0) - conform ipotezei prin care ecuaţia diferenţială ce formează problema Cauchy dă în orice punct (x,y) panta curbei!!!

Metodele Runge – Kutta

Metodele lui Euler implică necesitatea evaluării derivatelor de ordin superior ale funcţiei y(x) respectiv ale funcţiei f(x,y) care duc la dificultăţi în aproximarea numerică a derivatelor de ordin superior.

În schimb metodele de tip Runge – Kutta evită în totalitate utilizarea derivatelor de ordin superior ele folosind numai derivatele de ordin 1 ale funcţiei y(x), adică valorile funcţiei f(x,y).

Se calculează valorile funcţiei f(x,y) într-un număr de puncte intermediare ale intervalului [xi-1,xi] pentru determinarea lui yi cu o eroare minimă.

Cu alte cuvinte metodele Runge – Kutta de integrare numerică a unei ecuaţii diferenţiale, înlocuiesc calculul derivatelor funcţiei f(x,y) prin evaluări ale sale în diverse puncte.

Fie ecuaţia diferenţială ordinară cu condiţii iniţiale de forma:( )

( )⎩⎨⎧

==

00 yxyy,xf'y

hixx 0i ⋅+= ni ,...,2,1=Fie

o diviziune echidistantă a intervalului [x0, b]!!!

.

bxxx n =<<< ...10

Din raţiuni de simplificare a calculelor considerăm combinaţii liniare de valori ale funcţiei în anumite puncte ale intervalului [xi, xi-1], soluţia calculându-se cu o relaţie unipas de forma:

nn1100i1i ka...kakayy ⋅++⋅+⋅+=+

unde, din condiţia ca dezvoltarea în serie Taylor a membrului drept (în funcţie de h) să coincidă cu membrul drept al formulei lui Taylor de ordinul (n+1), avem şi formula de mai jos şi toţi coeficienţii după particularizări:

( )( )( )

( )1n1n,n11n00ninin

121020i2i2

010i1i1

ii0

k...kky,hxfhk...............................................................

kky,hxfhkky,hxfhk

y,xfhk

−− ⋅λ++⋅λ+⋅λ+⋅μ+⋅=

⋅λ+⋅λ+⋅μ+⋅=⋅λ+⋅μ+⋅=

⋅=nn1100i1i ka...kakayy ⋅++⋅+⋅+=+

Particularizând parametrul n determinăm diverse formule de tip Runge – Kutta:

1. n=0 (Runge – Kutta de ordin I), y0 - dat( )iii1i y,xfhyy ⋅+=+

2. n=1 (Runge – Kutta de ordin II), y0 - dat ( )10i1i kk21yy ++=+

( )ii0 y,xfhk ⋅= ( )0ii1 ky,hxfhk ++⋅=

( ) ( )( )[ ]iiiiiii1i y,xfhy,hxfy,xf2hyy ⋅++++=+

- formula modificată a lui Euler.

- formula lui Euler clasică

3. n=2 (Runge – Kutta de ordin III ), y0 – dat ( )210i1i kk4k61yy +++=+

( )ii0 y,xfhk ⋅= ⎟⎠

⎞⎜⎝

⎛ ++⋅=2

ky,

2hxfhk 0

ii1

( )01ii2 kk2y,hxfhk −++⋅=

( )3210i1i kk2k2k61yy ++++=+

( )ii0 y,xfhk ⋅= ⎟⎠

⎞⎜⎝

⎛ ++⋅=2

ky,

2hxfhk 0

ii1

⎟⎠

⎞⎜⎝

⎛ ++⋅=2

ky,

2hxfhk 1

ii2 ( )2ii3 ky,hxfhk ++⋅=

4. n=3 (Runge – Kutta de ordin IV), y0 – dat

Este foarte utilizată în aplicaţiile din domeniul electrotehnic, complicate şi pretenţioase din punct de vedere a preciziei!!!

Metode numerice indirecte ( metode multipas )

Soluţia în punctul xi se determină prin aceste metode multipas folosindu-se valorile calculate ale funcţiei în mai mulţi paşi anteriori. Dacă la metodele unipas funcţia necesită evaluări pentru un număr mare de valori ale variabilei independente la metodele multipas (metode cu paşi legaţi) nu este necesar calculul valorilor funcţiei f(x,y) în puncte intermediare suplimentare faţă de cele corespunzătoare pasului de discretizare (integrare) h. Pentru aceste metode este preferabil ca punctele luate în considerare pentru calculul soluţiilor sa fie echidistante.

Adams Adams – Bashforth Adams – Moulton

Predictor – corector Milne Predictor - corector Hamming

Fiind date valorile (x0,y0), (x1,y1),..., (xi,yi), metodele multipas folosesc aceste informaţii pentru calculul lui yi+1. Dezavantajul acestor metode este pornirea mai dificilă. Ele nu se autopornesc, adică la primul pas nu sunt disponibile (nu se cunosc) informaţiile din punctele anterioare necesare (adică primele valori ale soluţiei trebuie să fie calculate prin alte metode)!!!

Metoda dezvoltarii in serie Taylor

Fie ecuatia diferentiala ordinara cu conditia initiala Cauchy pusa:

y' x 0.1 y⋅+ y 0( ) 1 x ia valori pe [0;1]

f x y,( ) x 0.1 y⋅+:= \\ functia asociata membrului drept al ecuatiei explicite

\\ marginile intervalului de definitiesi numarul de puncte de calcula 0:= b 3:= N 100:= h

b a−

N:=

y0 1:= \\ conditia initiala a ecuatiei diferentiale de ordinul I

i 0 N..:= xi a h i⋅+:= \\ relatia de salt in nodurile de calcul

fx x y,( )xf x y,( )d

d:= fy x y,( )

yf x y,( )d

d:= \\ derivatele partiale ale functiei

yi 1+ yi h f xi yi,( ) h2

fx xi yi,( ) f xi yi,( ) fy xi yi,( )⋅+( )⋅+⎡⎢⎣

⎤⎥⎦

⋅+:= \\ formula de aproximare a functiei necunoscute;

yT 0 1 2 3 4 5 6 7 8 9

0 1 1.003 1.008 1.013 1.019 1.026 1.034 1.043 1.053 1.064=

0 0.5 1 1.5 2 2.5 30

5

10

y

x

Metoda lui Euler clasica (a liniilor poligonale)

Fie ecuatia diferentiala ordinara cu problema Cauchy pusa:

y' x 0.1 y⋅+ y 0( ) 1 x ia valori pe [0;1]

f x y,( ) x 0.1 y⋅+:= \\ functia asociata membrului drept al ecuatiei explicite

a 0:= b 1:= N 100:= hb a−

N:= \\ marginile intervalului de definitie

si numarul de puncte de calcul

y0 1:= \\ conditia initiala a ecuatiei diferentiale de ordinul I

i 0 N..:= xi a h i⋅+:= \\ relatia de salt in nodurile de calcul

yi 1+ yi h f xi yi,( )⋅+:= \\ formula de aproximare

yT 0 1 2 3 4 5 6 7 8

0 1 1.001 1.002 1.003 1.005 1.006 1.008 1.009 1.011=

\\ rezultate numerice

0 0.2 0.4 0.6 0.8 11

1.2

1.4

1.6

1.81.617

1

y

10 x

1. Sa se studieze conectarea unui circuit alcatuit dintr-o rezistenta r si o inductanta L la o sursa de tensiune continua U si scurtcircuitarea circuitului.

a) Conectarea circuitului la sursa continua

Legea a doua a lui Kirchhoff devine: U r i⋅ Ltid

d⋅+

Curentul este format din doua componente una stationara si una tranzitorie.

i ist itr+

Ltid

d⋅ U r i⋅−

diU r i⋅−

dtL

U i r⋅− A e

r−

Lt⋅

⋅

Tinand cont de conditiile initiale rezulta: t 0 i 0 U A

i I 1 e

t−

τ−

⎛⎜⎝

⎞⎟⎠⋅

unde τ se numeste constanta de timp a circuitelor si se masoara in secunde: τ

Lr

Componenta stationara a curentului este: istUr

Componenta tranzitorie a curentului devine: itrU−

re

r−

Lt⋅

⋅

Valori numerice:

L 0.001:= [H] r 0.1:= [Ω] U 10:= [V]

Constanta de timp a circuitului este: τLr

:= τ 0.01= [s]

Componenta stationara a curentului este: ist t( )Ur

:=

Componenta tranzitorie a curentului este: itr t( )U−

re

t−

τ⋅:=

Curentul total:

i t( ) ist t( ) itr t( )+:= t 0 0.0001, 0.08..:= tau t( )ist t( )

τt⋅:=

0 0.02 0.04 0.06 0.08

100

100

comp stationaracomp tranzitoriecurentul totalconst de timp

τ

Tensiunea sinusoidala la care este alimentat circuitul este de forma:u Um sin ω t⋅ ψ+( )⋅

Curentul este de forma: i ist itr+ ist A e

r−

Lt⋅

⋅+

0 0.01 0.02 0.03 0.04 0.05 0.06

50

50

comp stationaracomp tranzitoriecurentul total

42.082

30.331−

i st t ψ,( )

i tr t ψ,( )

i t ψ,( )

0.060

τ

t

2. Circuitul din figura de mai jos functioneaza in regim permanent cu intrerupatorul K deschis. Sa se determine variatia in timp a curentului din bobina, in urma scurtcircuitarii rezistorului R.1. Se precizeaza valorile numerice ale parametrilor circuitului si a sursei de alimentare.

=>

R 10:= Ω R1 30:= Ω L 0.1:= H --> parametrii circuitului electric

Semnalul furnizat de sursa de alimentare se va defini de tip ferastrau periodic cu posibilitatea de a modifica alura de crestere, dupa cum urmeaza:

0 2 4 6 8 10

25

50

75

100

125

150100.323

0

u t( )

100 t

A. Modificarea pozitiei intrerupatorului K conduce la aparitia regimului tranzitoriu in circuitul electric R-L. Modelul matematic, reprezentat de o ecuatie diferentiala, se obtine din aplicarea teoremei a doua a lui Kirchhoff pe ochiul de circuit. Conditia initiala a ecuatiei diferentiale se deduce din calculul intensitatii curentului electric in regimul permanent anterior aparitiei fenomenului tranzitoriu. Se modifica succesiv forma tensiunii de alimentare pentru a se vedea tipul de variatie al curentului electric prin circuit.

B.L

tid

d⋅ R i⋅+ u t( ) (ecuatia diferentiala de ordinul I, numita)

T0 2.5:= sec (perioada de timp la care se face comutatia)

i0u T0( )

R R1+:= (conditia initiala la perioada T.o ,deoarece regimul

de functionare nu este alternativ, ci doar periodic).

i0 1.438= A

0 C. Desi solutia analitica a ecuatiei diferentiale descrise anterior se cunoaste, mai intai se va solutiona problema printr-o metoda numerica de integrare aproximativa, numita Runge-Kutta de ordinul 4. Scopul acestei alternative este de a prezenta o varianta complet bazata pe calcule numerice la varianta tratarii analitice. Se scrie derivata curentului in functie de celelalte marimi din ecuatie; se ataseaza membrului drept al ecuatiei o functie de variabile timpul si curentul electric; se conditioneaza perioada de simulare, incepand de la comutarea intrerupatorului si un pas de discretizare a timpului de calcul; un set de formule numite Runge-Kutta intra intr-o formula recursiva de aproximare a formei de variatie a curentului electric.

tid

du t( )

LRL

i⋅− => f t i,( )u t( )

LRL

i⋅−:= (functia atasata)

f x y,( )u x( )

LRL

y⋅−:= (argumentele functiei s-au denumit altfel in sensul unei matematizari a problemei)

T0 2.5:= sec Tf 6:= sec h 0.01:= (pasul de discretizare a timpului)

MTf T0−

h:= M 350= (in dependenta directa cu numarul de puncte

de calcul se afla precizia metodei numerice de aproximare a valorilor de curent electric).

i 0 M..:= xi T0 h i⋅+:= --> punctele de timp de calcul din intervalul de simulare;

y0 i0:= --> initializarea conditiei impusa in ecuatia diferentiala;

k0 x y,( ) h f x y,( )⋅:= k1 x y,( ) h f xh2

+ yk0 x y,( )

2+,

⎛⎜⎝

⎞⎟⎠

⋅:=

k2 x y,( ) h f xh2

+ yk1 x y,( )

2+,

⎛⎜⎝

⎞⎟⎠

⋅:=

(setul de functii Runge-Kutta)

k3 x y,( ) h f x h+ y k2 x y,( )+,( )⋅:=

--> formula recursiva de aproximare a valorilor de curent electric in punctele de timp discretizate:

y i 1+ y i16

k0 xi y i,( ) 2 k1 xi y i,( )⋅+ 2 k2 xi y i,( )⋅+ k3 xi y i,( )+( )⋅+:=

2.5 3 3.5 4 4.5 5 5.5 6

2

4

6

8

10

y

i0

x T 0,

E. Problema de regim tranzitoriu expusa in acest caz nu prezinta complexitate, insa faptul ca la varianta analitica dedusa se contrapune o solutionare numerica, constituie un plus de cunoastere, de verificare si experimentare. Si mai mult, extinderea la circuite mai complexe in regim tranzitoriu nu aduce dificultati la alternativa numerica, pe cand cea analitica se complica.

1. Sa se rezolve (sa se determine solutia aproximativa) utilizand metoda lui

Euler ecuatia diferentiala y'1

1 x2+

y⋅ cu conditia initiala y 0( ) 1 pe intervalul

[0,1] si cu pasul h=0.02.

f x y,( )1

1 x2+

y⋅:= ti 0:= tf 1:= h 0.01:=

Ntf ti−

h:= N 100= n 1 N..:= r 0 5, N..:= y0 1:=

Punctul de pornire: x0 ti:= xn x0 n h⋅+:= yn yn 1− h f xn 1− yn 1−,( )⋅+:=

Ultimul punct calculat: xN 1= yN 2.192=

r05

10

15

20

25

30

35

40

45

50

55

60

65

70

75

= xr

00.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0.55

0.6

0.65

0.7

0.75

= yr

11.051

1.104

1.16

1.217

1.276

1.337

1.399

1.461

1.525

1.588

1.652

1.715

1.778

1.84

1.901

=

0 0.5 11

1.5

2

2.52.192

1.01

y n

10.01 xn

Rezolvam aceeasi ecuatie diferentiala cu ajutorul functiei de rezolvare rkfixed predefinita in Mathcad, care utilizeaza pentru rezolvare metoda Runge-Kutta IV.

f x y,( )1

1 x2+

y⋅:= ti 0:= tf 1:=

N 1000:= i 0 100, N..:= val 0 1:=

D x Y,( ) f x Y0,( ):= S rkfixed val ti, tf, N, D,( ):= X S 0⟨ ⟩:= Y S 1⟨ ⟩

:=

Xi

00.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

= Yi

11.105

1.218

1.338

1.463

1.59

1.717

1.842

1.964

2.081

2.193

=Y N 2.193= i0

100

200

300

400

500

600

700

800

900

1·10 3

=

0 0.5 11

1.5

2

2.5

Yi

Xi

2. Sa se rezolve ecuatia diferentiala y' 1 x− 4 y⋅+= cu conditia initiala y(0)=1 si cu pasul h=0.01 pe intervalul [0,1], utilizand metoda Runge - Kutta de ordinul IV.

f x y,( ) 1 x− 4 y⋅+:= ti 0:= tf 1:= h 0.01:=

Ntf ti−

h:= N 100= n 1 N..:=

Conditia initiala: y0 1:= x0 ti:= xn x0 n h⋅+:=

Notam:

K1 x y,( ) f x y,( ):= K2 x y,( ) f xh

2+ y

h

2K1 x y,( )⋅+,⎛⎜

⎝⎞⎟⎠

:=

K3 x y,( ) f xh

2+ y

h

2K2 x y,( )⋅+,⎛⎜

⎝⎞⎟⎠

:= K4 x y,( ) f x h+ y h K3 x y,( )⋅+,( ):=

y n yn 1−

h

6K1 x

n 1−y

n 1−,( ) 2 K2 x

n 1−y

n 1−,( )⋅+

2 K3 xn 1−

yn 1−

,( )⋅ K4 xn 1−

yn 1−

,( )++

...⎛⎜⎜⎝

⎞⎟⎟⎠

⋅+:=yN 64.898=

r 0 10, N..:=r

010

20

30

40

50

60

70

80

90

100

= xr

00.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

= y r

11.609

2.505

3.83

5.794

8.712

13.053

19.516

29.145

43.498

64.898

=

0 0.5 10

50

100

yr

xr

Rezolvam aceeasi ecuatie diferentiala cu ajutorul functiei de rezolvare rkfixed predefinita in Mathcad, care utilizeaza pentru rezolvare metoda Runge-Kutta IV.

ti 0:= tf 1:= N 1000:= i 0 100, N..:= val0 1:=

D x Y,( ) f x Y0,( ):= S rkfixed val ti, tf, N, D,( ):= X S0⟨ ⟩

:= Y S1⟨ ⟩

:=

YN 64.898= Tabelele valorilor calculate:

0 0.5 10

20

40

60

80

Yi

X i

i

0100

200300

400

500

600

700

800

900

1·10 3

= X i

00.1

0.20.3

0.4

0.5

0.6

0.7

0.8

0.9

1

= Yi

11.609

2.5053.83

5.794

8.712

13.053

19.516

29.145

43.498

64.898

=

7. Scriem sistemul urmat in partea dreapta de conditiile init

xy0 x( )d

d5− y0 x( )⋅ y1 x( )+ y2 x( )+ y0 0( ) 1

ty1 x( )d

d3 y0 x( )⋅ 2 y0 x( )⋅ y2 x( )⋅− y1 0( ) 0

ty2 x( )d

d2 y0 x( )⋅ y1 x( )⋅ 6 y2 x( )⋅− y2 0( ) 1−

D x Y,( )

5− Y0⋅ Y1+ Y2+

3 Y0⋅ 2 Y0⋅ Y2⋅−

2 Y0⋅ Y1⋅ 6 Y2⋅−

⎛⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎠

:= Y0

1

0

1−

⎛⎜⎜⎝

⎞⎟⎟⎠

:=

x0 0:= Valoarea initiala a variabilei independentex1 10:= Valoarea finala a variabilei independente

N 500:= Numarul de valori ale solutiei in [x0, x1]Matricea solutiei: S Rkadapt Y0 x0, x1, N, D,( ):=

Valorile variabilei independente x S 0⟨ ⟩:=

ST

0 1 2 3 4 5 6 7 8 9 10

01

2

3

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.21 0.888 0.79 0.705 0.631 0.566 0.51 0.461 0.418 0.381 0.349

0 0.092 0.17 0.237 0.294 0.344 0.386 0.424 0.457 0.486 0.512

-1 -0.885 -0.781 -0.687 -0.603 -0.527 -0.46 -0.401 -0.348 -0.302 -0.261

=

xT 0 1 2 3 4 5 6 7 8 9

0 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18=

y1T 0 1 2 3 4 5 6 7 8

0 0 0.092 0.17 0.237 0.294 0.344 0.386 0.424 0.457=

y2T 0 1 2 3 4 5 6 70 -1 -0.885 -0.781 -0.687 -0.603 -0.527 -0.46 -0.401

=

y0T 0 1 2 3 4 5 6 7 8

0 1 0.888 0.79 0.705 0.631 0.566 0.51 0.461 0.418=

8. Metoda Euler pentru rezolvarea ecuatiilor diferentiale

n 10:= (numarul de subintervale ale intervalului [0,1]) h1n

:= k 0 n..:=

y0 1.5:= x0 0:= (datele initiale ale metodei)

xk x0 k h⋅+:= (nodurile utilizate in metoda Euler)

f x y,( ) 1 x y2⋅+:= (functia din membrul drept al ecuatiei diferentiale)

yk 1+ yk h f xk yk,( )⋅+:= formula iterativa - calculeaza solutia pe noduri (solutia discreta!)

yT 0 1 2 3 4 5 6 7 8 9

0 1.5 1.6 1.726 1.885 2.092 2.367 2.747 3.3 4.162 5.647=

Reprezentarea grafica a solutiei determinata cu metoda lui Euler

0 0.5 10

5

10

yk

xk

O solutie continua (aproximativa) se determina, de exemplu, prin interpolare Lagrange:

L t( )

k 0

n

i

if i k 1,t xi−

xk xi−,

⎛⎜⎝

⎞⎟⎠∏

=

⎛⎜⎜⎝

⎞⎟⎟⎠

yk⋅∑:= L xk( )1.51.6

1.726

1.885

2.092

2.367

2.747

= yk1.51.6

1.726

1.885

2.092

2.367

2.747

=

0 0.5 10

5

108.618

1.5

L t( )

y

10 t x,