Dan.Micu@et.utcluj.ro

Curs 2. Metode numerice de rezolvare a ecuaţiilor

neliniare algebrice şi transcendente

Aplicaţii în Ingineria Electrică

METODE NUMERICE

Ecuaţiile neliniare constituie una din cele mai frecvente aplicaţii decalcul numeric care apare în cadrul activităţilor de proiectare diningineria electrică:

Rezolvarea numerică a ecuaţiei de stare a conductorului LEAAnaliza stabilităţii unui generator sincron prin localizarea rădăcinilor ecuaţiei caracteristice

Problema…Expresia ecuaţiei este foarte complicată sau valorile coeficienţilor nu

se cunosc exact (rezultaţi din determinări experimentale sau au fost calculaţi pe baza unor ipoteze simplificatoare)!

Concluzia…Nu se pune problema soluţionării exacte a ecuaţiilor cu metode directe (nr. finit de

paşi – cunoscut apriori)Se utilizează metode numerice aproximative – metode iterative cu convergenţă

teoretic infinită dar practic finită (prin estimarea permanenta a gradului de precizie a determinării soluţiei sau soluţiilor)

Soluţiile ecuaţiilor neliniare se obţin aşadar ca limite ale unor şiruri convergente.

Consideraţii teoretice

Se consideră o ecuaţie de formă generală:( ) 0xf =

Pentru majoritatea aplicaţiilor din ingineria electrică – domeniul de definiţie este un interval I:

RI:f →

Teorema lui Rolle-de localizare a rădăcinilor:

( ) ( ) ( )b,ax0bfaf 0 ∈∃⇒<⋅

x0-rădăcina ecuaţiei

Observaţii

Dacă f(x) este un polinom atunci ecuaţia se numeşte ecuaţie algebrică; în caz contrar se numeşte ecuaţie transcendentă

Se numeşte: - rădăcina funcţiei (ecuaţiei) – numai la ecuaţii algebrice

- zeroul (soluţia) funcţiei - la ecuaţii transcendente

RezumatSe va demonstra pas cu pas modul în care o aplicaţie reală

(provenită din ingineria electrică) se modelează matematic şi apoi

se rezolvă numeric cu ajutorul metodelor numerice;

Intuitiv, vor fi experimentate fazele de soluţionare numerică a

ecuaţiilor algebrice (polinomiale), stabilite ca model matematic al

aplicatiei reale;

De la desfăşurarea particularizată a metodei numerice de

rezolvare, se trece la descrierea ei generală, fiind subliniate limitele

de aplicabilitate, avantaje şi dezavantaje;

Societatea de transport a energiei electrice, S.C TRANSELECTRICA S.A, şi-a stabilit ca obiectiv de investiţii pe anul 2009 construirea unei linii electrice de interconexiune între două staţii electrice. Pe distanţa celor două staţii, datorită diferenţelor de teren, amplasamentul se împarte în două zone, delimitate printr-o linie WE, aşa încât costul execuţiei liniei pe fiecare zonă se caracterizează prin indicii C1 şi C2.

1. Metoda bisecţiei (a înjumătăţirii intervalului)

Se pune problema stabilirii traseului optim al liniei astfel încât costul de execuţie să fie minim.

Pentru soluţionarea modelului matematic care provine din aceastaaplicatie este necesar apelul la o metodă numerică.

Aplicaţie practică

Formularea Problemei (P) – date cunoscute (date); necunoscute (soluţii); lege de legătură (date-soluţii)

Descrierea Problemei (P) – Model Matematic (M(P))

Aproximare M(P) – printr-o Metodă Numerică (MN(P))

Implementare algoritm în MathCad (Matlab, Mathematica etc.)

Dezvoltare algoritmPentru MN(P)

Specificareaproblemei reale

(Traseu optim LEA)

Construirea unui model fizic(geometria traseului)

Confruntarea cu realitatea

Formularea problemeimatematice

(ecuaţie algebrică)

Interpretarea soluţiei (cost minim)

Rezolvarea problemei matematice

inducere

deducere

Metode numerice(Bisecţia)

Modificare model

Fig. 1 Reprezentarea geometrică de amplasament a liniei electrice

Ideea de bază se reduce la localizarea unui punct P pe frontiera WE prin care linia electrică să traverseze limita de separaţie dintre cele două zone. Având la dispoziţie indicii de cost, se poate stabili o relaţie de egalitate între sinusurile unghiurilor formate de direcţiile traseelor de linie din cele două zone şi perpendiculara dusă prin punctul P la zona de delimitare:

2211 sinsin

Modelul matematic pe baza datelor problemei tehnice de soluţionat:

θθ ⋅=⋅ CC( )( )22

22

222

22

1 xLbxLC

xaxC

−+−

⋅=+

⋅

( )[ ] ( ) ( )22222

22221 xaxLCxLbxC +⋅−⋅=−+⋅⋅ - ecuaţie algebrică

Exprimarea devine mai sugestivă prin înlocuirea datelor cunoscute:

a=3 km; b=1 km; WE=L=12 km; C1=87000 RON/km; C2=93000 RON/km.

0112091041868184225792259201080 234 =−⋅+⋅−⋅+⋅− xxxx

Necunoscuta în această ecuaţie polinomială (algebrică) este distanţa x de la marginea W la punctul P (capetele se consideră ştiute xW şi xE).

Dacă membrul stâng al acestei ecuaţii se exprimă ca o funcţie f(x), evident derivabilă pe intervalul [xW;xE], se permite efectuarea următoarelor testări:

Coeficienţi determinaţi experimental şi prin simplificări în model!!!

112091041868184225792259201080)( 234 −⋅+⋅−⋅+⋅−= xxxxxf

Intre cele două valori limită pe care poate să le ia necunoscuta x trebuie să existe o valoare care să anuleze funcţia (Rolle)!!!

0)x(f)x(f EW <⋅Dacă:

Reprezentarea grafică a polinomului pe intervalul definit arată clar îndeplinirea condiţiei anterioare:

Fig. 2 Variaţia polinomului şi partiţionarea intervalului de definire

Dar mai mult decât atât, în vederea localizării cât mai exacte a soluţiei, se sugerează ideea împărţirii intervalului prin înjumătăţire succesivă:

( )EW xxc +⋅=21

0)()( <⋅ cfxf W

0)()( =⋅ cfxf W

0)()( >⋅ cfxf W

Dacă soluţia se află în intervalul [xW ; c];

soluţia este chiar valoarea c;

soluţia se află în cealaltă jumătate de interval.

Observaţie: reprezentarea grafică a funcţiei oferă ca evidentă validitatea celei de-a treia ipoteze.

( )Excd +⋅=21

0)()( <⋅ Exfdf

Se înjumătăţeşte intervalul [c; xE] prin aceeaşi formulă de mediere aritmetică:

Variaţia grafică a funcţiei polinomiale arată clar verificarea inegalităţii !!!

Indică localizarea soluţiei în intervalul delimitat de punctele în care se face această evaluare!!!

Fig. 3 A doua înjumătăţire a intervalului

Continuând cu înjumătăţirea intervalului, se observă, în figura 4, restrângerea domeniului în care se află soluţia, apropierea tot mai certă de valoarea care anulează funcţia polinomială!

Fig. 4 Continuarea partiţionării intervalului; apropierea de soluţie

Procesul de restrângere treptată a intervalului de definire se poate derula până când evaluarea funcţiei în variabila de înjumătăţire (c, d, ...) devine mai mică decât o valoare impusă, ori efectiv se anulează. În oricare din aceste situaţii, se consideră determinată soluţia realizabilă a polinomului în limita unei precizii, impuse apriori.

Modalitatea prin care s-a soluţionat aplicaţia propusă nu reprezintă altceva decât o metodă numerică, numită metoda bisecţiei (metoda înjumătăţirii intervalului).

Pornind de la cele expuse, se va căuta descoperirea caracterului de generalitate al acestei metode.

Ideea se transpune generalizat după următorul algoritm:

];[];[ 00 baba =

0=k

( )kkkk abac −⋅+=+ 21

1

1111 ;0)()( ++++ ==⇒<⋅ kkkkkk cbaaafcf

kkkkkk bbcabfcf ==⇒<⋅ ++++ 1111 ;0)()(

1+= kk

Pasul 1: Se iniţializează limitele intervalului de căutare în care se caută soluţia cu Rolle

Pasul 2: start iterativ

Pasul 3: la un pas oarecare al procesului de calcul se determină noua valoare a soluţiei

Pasul 4: La acelaşi pas se calculează f (ck+1), f(ak) rezultând noile limite ale intervalului de căutare:

Pasul 5: Dacă:

Pasul 6: Incrementează şi reia Pasul 3.

Oricât de mult s-ar restrânge intervalul în jurul soluţiei, există posibilitatea ca valoarea determinată considerată drept soluţie să nu fie cea adevărată.

Mai mult, procesul iterativ de înjumătăţire nu poate continua la infinit, ci trebuie oprit după un anumit număr de partiţionări. Se va determina acest număr prin stabilirea unei erori limită între valoarea determinată ca şi soluţie şi valoarea adevărată.

Fie ε o limită de eroare absolută impusă.

( ) ( ) ( )00

3

11

2

2233 21

21

21 abababab −⋅⎟

⎠⎞

⎜⎝⎛=−⋅⎟

⎠⎞

⎜⎝⎛=−⋅=−

sau prin generalizare

( ) ( )0011 21...

21 ababab

n

nnnn −⋅⎟⎠⎞

⎜⎝⎛==−⋅=− −−

unde n este numărul până la care a ajuns iteraţia k.

Se ştie şi că eroarea cunoscută verifică inegalitatea: εα ≤− nxα – valoarea adevărată a soluţiei xn - valoarea aproximată la iteraţia n.

Însă diferenţa dintre valoarea adevărată şi valoarea de la a n-a iteraţie nu poate fi mai mare decât diferenţa dintre valorile capetelor intervalului curent restrâns:

( ) ( ) ( )abababxnn

nnn −⋅⎟⎠⎞

⎜⎝⎛=−⋅⎟

⎠⎞

⎜⎝⎛=−≤−

21

21

00α ( ) ε≤−⋅⎟⎠⎞

⎜⎝⎛ ab

n

21 ( ) ( )εloglog

21log ≤−+⎟⎠⎞

⎜⎝⎛⋅ abn

( ) ( )( )2log

loglog ε−−≥

abn

Oferă o măsură destul de precisă a numărului minim de iteraţii necesare pentru calcularea unei soluţii aproximative care să se încadreze într-o limită de eroare impusaapriori!!!

⇒ ⇒

Este posibilă adăugarea unei condiţii de oprire a procesului când se atinge o anumită precizie, adică după ce se efectuează un număr de n iteraţii.

Concluzii asupra metodei bisectiei

Metoda bisecţiei converge spre soluţie indiferent cât de departe estepunctul de pornire, fapt pentru care se consideră o metodă globalăde determinare a soluţiilor;

Convergenţa spre soluţie este lentă, trebuie executat un numărmare de împărţiri ale intervalului pentru a ajunge la o precizie satisfăcătoare;

Metoda bisecţiei nu poate fi utilizată pentru determinarea soluţiilorunei funcţii care este tangentă la axa Ox, fără să o străpungă, fiindcă nu se verifică chiar condiţia de start.

0)()( <⋅ bfaf

Solutionarea aplicatiei in utilitarul Mathcad cu un algoritm al metodei bisectiei

Fie ecuatia polinomiala de gradul IV:

1080− x4⋅ 25920x3

⋅+ 225792x2⋅− 1868184x⋅+ 11209104− 0

Atasam acestei ecuatii functia corespondenta :

f x( ) 1080− x4⋅ 25920 x3

⋅+ 225792 x2⋅− 1868184 x⋅+ 11209104−:=

Reprezentam grafic functia pe domeniul de definitie:a 0:= b 12:= x a b..:=

0 2 4 6 8 10 12

2 .107

1 .107

1 .107

f x( )

x Fixam o precizie de calcul, adica o eroare absoluta cu care se va determina radacina de interes a ecuatiei, fata de valoarea adevarata: ε 10 6−:=

Executam desfasurat primii doi pasi ai metodei bisectiei:se initializeaza marginile intervalului de pornire a 0 a:= si b 0 b:= •

se testeaza conditia de existenta a solutiei in acest interval •f a0( ) f b0( )⋅ 0< 1=

se realizeaza prima injumatatire a intervalului•

apoi fixarea noilor margini ale intervalului restrans si testarea conditiei de existenta a solutiei:•

a2 c2:= b2 b1:= f a2( ) f b2( )⋅ 0< 1=

urmeaza a doua injumatatire a intervalului•

c2 a1

b1 a1−

2+:= c2 9= \\ a doua aproximare a solutiei

se testeaza conditia de existenta a solutiei f a1( ) f b1( )⋅ 0< 1=•

c1 a0

b0 a0−

2+:= c1 6= \\ prima aproximare a solutiei

se trece la noile margini ale intervalului, restrans, in care se afla solutia•

a1 c1:= b1 b0:=

Conform demonstratiei referitoare la numarul de iteratii necesare pentru realizarea unei precizii impuse, se evalueaza efortul minim de calcul:

n roundlog b0 a0−( ) log ε( )−

log 2( )

⎛⎜⎝

⎞⎟⎠

:= n 24=

In continuare, pentru obtinerea solutiei cautate se propune un algoritm de calcul automat, corespunzator celui dat in partea teoretica sub forma de pseudolimbaj:

a 0:= b 12:=

sol

a0 a←

b0 b←

ck 1+ ak

bk ak−

2+←

ak 1+ if f ak( ) f ck 1+( )⋅ 0< ak, ck 1+,( )←

bk 1+ if f ck 1+( ) f bk( )⋅ 0< bk, ck 1+,( )←

a b∧

k 0 n..∈for

c

:= \\ se apeleaza ciclul "for"\\ se initializeaza marginile; primului interval in care se afla solutia;

\\ formula de aproximare a solutiei;\\ modificarea la fiecare iteratie a capetelor intervalului aflat in restrangere;\\ memorarea in vectori a acestor margini;\\ algoritmul returneaza vectorul aproximatiilor solutiei, care sunt de fapt injumatatirile intervalelor.

Afisam aproximatiile realizate prin metoda bisectiei pentru zerourile functiei polinomiale:

soln 9.983287=

solT 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

0 0 6 9 10.5 9.75 10.125 9.938 10.031 9.984 9.961 9.973 9.979 9.981 9.983 9.984 9.983 9.983=

0 5 10 15 20

2468

1012

solk

solk

k

Reprezentam grafic procesul de convergenta a aproximatiilor catre valoarea considerata optima a solutiei, din punct dvedere al preciziei:

k 0 n..:=

Observatii: 1. procesul de convergenta solicita un efort de calcul considerabil, ceea ce poate fi un dezavantaj al metodei; 2. in practica, pentru problema propusa avem nevoie doar de 4 zecimale in rezultatul calculat; 3. algoritmul de mai sus poate fi imbunatatit prin impunerea unei conditii de stopare cand se atinge o anumita precizie

f t( ) g02 1 2 sin t( )⋅+ sin t( )2

+( )⋅ 2 g0⋅ v⋅ t⋅ 1 sin t( )+( )⋅− v2 t2⋅+⎡⎣

⎤⎦ iD⋅ ε U⋅ v⋅ S⋅−:=

g02 1 2 sin t( )⋅+ sin t( )2

+( )⋅ 2 g0⋅ v⋅ t⋅ 1 sin t( )+( )⋅− v2 t2⋅+⎡⎣

⎤⎦ iD⋅ ε U⋅ v⋅ S⋅− 0=

C. Cerinta problemei impune aducerea ecuatiei transcendente de la forma redata ca model matematic, corespunzator fenomenului electric din condensator, la o forma convenabila exprimarii ca ecuatie careia sa i se poata atasat o functie:

iD εU

g0 1 sin t( )+( )⋅ v t⋅−⎡⎣ ⎤⎦2

⋅ v⋅ S⋅=<=>iD εU

g t( ) v t⋅−( )2⋅ v⋅ S⋅=

iD S→

JD→⎯⌠

⎮⎮⌡

d==>JD→⎯

tE→

ε⋅∂

∂==>E

→ Ug t( ) v t⋅−

u→⋅=

B. Din expresia intensitatii campului electric, iar apoi din variatia inductiei electrice in timp se exprima analitic formula intensitatii curentului electric de deplasare intre armaturi:

A. Modificarea raportului metric dintre armaturi conduce la aparitia unui curent de deplasare dependent de legile impuse miscarii celor doua armaturi.

Solutie

--> functia dupa care se modifica spatiul dintre armaturig t( ) g0 sin t( ) 1+( )⋅:=

ε ε0 εr⋅:=Fm

ε01

4 π⋅ 9⋅ 109⋅

:=ms

v 50:=mg0 0.1:=

AiD 3 10 6−⋅:=εr 2.5:=VU 10000:=m2S 10 3−

:=S 10 cm2⋅:=

3. Un condensator plan are fiecare armatura de suprafata S; distanta initiala dintre armaturi este g0. Condensatorul este cufundat intr-o baie de ulei si este

conectat la diferenta de potential U. Intr-un moment, dat se deplaseaza o armatura catre cealalta cu viteza relativa v; in acelasi timp, cea de-a doua armatura se departeaza fata de pozitia initiala astfel incat distanta initiala dintre armaturi variaza dupa o functie data g(t). Se cere sa se calculeze timpul necesar ca valoarea curentului de deplasare dintre armaturi sa fie de 3 μA.

sol 0.000787= s Asadar, timpul la care curentul de deplasare ajunge la valoarea prescrisa este:

μs 10 6− s⋅:= tprescris 787 μs⋅:=

D. In utilitarul de calcul Mathcad exista o modalitate de rezolvare predefinita pentru ecuatia obtinuta, cu solutionare interna de tipul metodei expuse mai sus. In continuare, se va afisa solutia pe calea predefinita:

t 1 10 3−⋅:= s --> timpul aproximat

t'prescris root f t( ) t,( ):= --> functia predefinita de solutionare a ecuatiei neliniare

t'prescris 0.000765= s t'prescris 765 μs⋅:=

E. Observatie: Timpul calculat prin metoda predefinita difera de cel calculat prin metoda explicitata a bisectiei. Valoarea de referinta se considera cea data de functia predefinita. In aceasta situatie, rolul metodei explicitate devine demonstrativ. Totusi, eroarea care intervine, desi problema este conceputa cu implicatii teoretice, se situeaza la valori reduse.

Eroare t'prescris tprescris−:= Eroare 22μs=

O solutie a ecuatiei neliniare de mai sus, in cazul in care se defineste un interval de timp probabil pentru functia continua atasata se gaseste prin impartirea repetata a intervalului in doua parti egale, asa incat se formeaza un sir de iterare. Metoda numerica se numeste a injumatatirii intervalului, avand ca si corespondent geometric bisectia. Conditiile de aplicare a metodei determina in acest caz, prin incercari, si fixarea domeniului de timp la care sa se restranga definirea functiei.

ti 0:= s --> timpul initial tf 2 10 3−⋅:= s --> timpul final

f ti( ) f tf( )⋅ 0< 1= --> verificarea conditiei de existenta a solutiei pe intervalul ales

Se alege un numar de iteratii pentru sirul care se formeaza prin injumatatiri repetate; conditionat, rezulta un pas de divizare si un set de puncte de divizare:

N 103:= i 0 N..:= h

tf ti−

N:= h 2 10 6−

×= tditi h i⋅+:= tinvi

tf h i⋅−:=

Algoritmul de calcul al valorilor sirului pana la obtinerea solutiei este redat in liniile urmatoare de program:

sol

A12

tditinvk

+⎛⎝

⎞⎠

⋅← f tdi⎛⎝

⎞⎠

f tinvk⎛⎝

⎞⎠

⋅ 0<if

k 0 N i−..∈for

i 0 N..∈for

A

:=

S.C ARMĂTURA S.A, firmă cu profil electromecanic, a primit o comandă de prelucrare a unor plăci metalice utilizate în construcţia releelor de telefonie mobilă. Operaţia principală executată asupra acestor piese constă în vopsirea prin pulverizare în câmp electric, procedeu numit vopsire electrostatică.

Problema care se pune în cadrul acestei aplicaţii se referă la găsirea unei legături între dimensiunile de vopsit ale unei plăci şi mărimile electrice prin care se poate ajusta procesul. Astfel, prin reglajul acestor mărimi, tensiune electrică, respectiv curent electric care parcurge plăcile, devine posibil controlul automat al suprafeţei de vopsit pentru fiecare plăcuţă.

Cunoscând caracteristicile şi dimensiunile instalaţiei de vopsire electrostatică, aplicaţia se poate modela simplu printr-un condensator plan în care se introduc simultan două piese de vopsit. În primul rând trebuie determinată distanţa maximă de pătrundere a plăcilor în interiorul condensatorului plan, distanţă care limitează direct suprafaţa ce urmează a fi acoperită cu vopsea a acestora.

2. Metoda lui Newton (metoda tangentei)Aplicaţie practică

Modelul prin care se poate caracteriza problema:

Fig. 6 Modelarea pieselor şi a instalaţiei de vopsire electrostatică

Dimensiunile geometrice ale sistemului condensator plan şi piese de vopsit:

L = 0.65 m ; l = 0.3 m ; l0 = 0.25 m ; d= 1 m ; g = 0.002 m ; y = 0.15 m ;

Mărimile electrice de alimentare:

U = 20000 V ; i = 3 A.

Rămâne să se determine ecuaţia de echilibru a forţei electrice cu cea electrodinamică, iar din această relaţie să se găsească mărimea variabilei x de pătrundere a pieselor în interiorul condensatorului.

Fig. 8 Piesele de vopsit – vedere de sus

Forţa de respingere dintre cele două piese se exprimă integral cu formula:

( ) ∫ ∫+−

−

+−

−+

⋅⋅⋅⋅⋅

=0 02

2

2

2

2121

20

20 1

2

l l

l

lxL

xL

xL

xLm dxdx

xxi

xFπ

μ

( ) ( ) 0=− xFxF me

( ) ( ) ( )( ) ( ) ( ) 0896.1295.234013ln28.17616.5

209ln55.1556.344023ln)28.17936.9(=−⋅+⋅−⋅⋅−+

+⋅−⋅−⋅+⋅−⋅⋅−xxx

xxxx

Iar din ecuaţia de echilibru a forţelor

şi după înlocuirile numerice se deduce o ecuaţie transcendentă:

În această ecuaţie necunoscuta este distanţa maximă de pătrundere a pieselor în interiorul condensatorului, în condiţiile în care sunt fixate mărimile electrice de alimentare.

Pentru a determina valoarea variabilei x care verifică ecuaţia neliniară de mai sus, ne vom folosi de următoarele transformări:-se notează membrul stâng al ecuaţiei cu o funcţie

( ) ( ) ( ) ( )( ) ( ) ( )896.1295.234013ln28.17616.5

209ln55.1556.344023ln)28.17936.9(−⋅+⋅−⋅⋅−+

+⋅−⋅−⋅+⋅−⋅⋅−=xxx

xxxxxf

-se dezvoltă în serie Taylor funcţia f(x) în jurul unui punct x0 reţinând doar primii doi termeni

( ) ( ) ( ) ( )00'

0 xxxfxfxf −⋅+=

-considerând pe x0 ca o aproximaţie iniţială a soluţiei ecuaţiei de la care s-a pornit, dacă în expresia dezvoltării Taylor de mai sus se înlocuieşte variabila x cu o nouă aproximare a soluţiei, x1, pentru care se presupune că f(x) se anulează, atunci rescriem:

( ) ( ) ( ) ( )( )0

'0

01010'

00xfxf

xxxxxfxf −=⇒−⋅+=

-în acest fel aproximaţia x1 devine calculabilă în raport cu prima aproximaţie, şi mai departe pentru a afla o aproximaţie cu o precizie sporită, efectuăm succesiv calculele

( )( )

( )( )n

nnn xf

xfxx

xfxfxx '1

1'

112 .... −=⇒⇒−= +

- ajungem astfel la o aproximaţie xn+1 a soluţiei, pe care în funcţie de numărul de iteraţii parcurse o vom adopta ca fiind soluţia căutată a problemei.

Numeric, pornind de la o aproximaţie iniţială x0 = 0.2 m, se ajunge după n = 9 iteraţii la x9 = 0.3 m.

Modalitatea prin care s-a soluţionat ecuaţia neliniară este aplicarea metodei numerice a lui Newton.

Caracterul general al metodei lui Newton

( ) ( ) ( ) ( )

Fie o ecuaţie de forma f(x)=0, cu variabila x din [a,b], iar funcţia continuă şi de două ori derivabilă pe intervalul dat. Dezvoltarea în serie Taylor în jurul unui punct xn, în cazul în care se reţin doar primii doi termeni ai dezvoltării şi restul, este:

( ) ( ) [ ]nnnnnnn xxcufxxxfxxxfxf ;!2

1 ''2' ∈⋅−⋅+⋅−+= ξξ

Înlocuirea în expresia de mai sus, a unei aproximaţii xn+1 în locul lui x, şi succesivă lui xn, pentru care presupunem că se anulează f(x) conduce la relaţiile:

( )( )n

nnn xf

xfxx '1 −=+

S-a generat astfel o formulă de calcul a soluţiilor ecuaţiilor, în care apare o eroare de metodă datorată neglijării restului seriei Taylor şi care se bazează pe un calcul recursiv.

( ) ( )( )

( ) ( )( )n

'n

''2

nn

'n

nxf

fxx21

xfxfxx0xf ξ

⋅−⋅−−=⇒=

neglijarea restului

1nxx +=⇒

Dacă se admite că reprezentarea grafică a funcţiei arată ca în figura 9, iar derivata lui f(x) se reprezintă ca o dreaptă care trece succesiv prin punctele nxxxx ...,,,, 210ALGORITM:1. Se iniţializează soluţia cu x0 (b=x0)2. La un pas oarecare k al procesului iterativ se calculează

f(xk), f ’(xk) rezultând noua valoare a soluţiei aproximative xk+1

3. Calculul se consideră terminat dacă:

Metoda lui Newton poate fi interpretată geometric: trasare repetată a tangentelor la funcţie!

Fig. 9 Interpretare geometrică

Procesul continuă până când se ajunge, în limita unei precizii, la soluţia considerată optimă.

Din start, în deducerea metodei lui Newton se introduce o eroare de aproximare prin neglijarea restului seriei Taylor.

Se va investiga mărimea ponderii pe care această eroare o are, modul prin care se caracterizează convergenţa procesului iterativ de aproximare şi efortul computaţional aferent acestui proces!!!

ε≤−+ k1k xx

Fie funcţia ( ) .,0),(2 solutieIptfsiRIICf ∈=⊂⊂ αα

Pentru un xn dat se defineşte:( )( )n

nnn xf

xfxx '1 −=+

Atunci, conform teoremei valorii medii, există un punct ξn între α şi x astfel încât din seria Taylor:

Interpretarea relatiei: eroarea la fiecare pas iterativ depinde de pătratul erorii la pasul anterior.

Convergenţa metodei este condiţionată: depinde de soluţia iniţială, recomandându-se pentru reducerea numărului de iteraţii, satisfacerea condiţiei:

Formula de oprire a procesului iterativ, confirmată numeric şi grafic de faptul că după un anumit număr de iteraţii, soluţia calculată trebuie să se stabilizeze:

( ) 51ε≤−+ −nnn xxxf

( ) ( ) 0x''fxf 00 >⋅

} ( )( )

} ( )( )n

'n

''2

n

x

n'

nn xf

fxx21

xfxfxx

1n

ξ⋅⎟

⎟

⎠

⎞

⎜⎜

⎝

⎛−⋅−−=

αα+48476

( ) ( ) ( )( )n

nnn xf

fxx '

''2

1 21 ξαα ⋅−⋅−=− +

Demonstrarea algoritmului Metodei lui Newton

Presupunem că x0, x1,... sunt aproximaţiile succesive

obţinute din procesul iterativ:

Se demonstrează că şirul este convergent şi limita şirului este soluţia ecuaţiei f(x)=0!!!

Presupunem ( ) ( ) 0x''f;0bf];b,a[x >>∈

Se alege x0=b şi deci: α>0x α - soluţia exactă

La pasul k presupunem: ( ) ( ) Nk,0x''f;0xf;x kkk ∈∀>>α>

Dacă: ( )kk xx −α+=α

( ) ( )( ) ( ) ( ) ( ) ( ) ( )444 3444 21kk x;stRe

k

2k

kk

k

Taylor

kk ''f!2xx'f

!1xxfxxf0f

<ξ<α

ξ⋅−α

+⋅−α

+=−α+==α

( ) ( ) ( ) ( ) ( )321

0k

2k

kkk ''f!2xx'fxxf

>

ξ⋅−α

−=⋅−α+

( )( )k

'k

k1k xfxfxx −=+

( ) ( ) ( ) ( ) ( ) ⇒<ξ⋅−α

−=⋅−α+ 0''f!2xx'fxxf k

2k

kkk( )( ) ⇒−<−α

k

kk x'f

xfx( )( )

48476 1kx

k

kk x'f

xfx

+

−<α

⇒α>>⇒⎭⎬⎫

><α

++

+1kk

1kk

1k xxxx

x Şirul este monoton descrescător şi mărginit de αdeci şirul este convergent (limită finită şi unică)

Presupunem că valoarea acestei limite este: kk

xlim∞→

=α

( )( )

( )( ) ( ) 0f

'ff

xfxfxlimxlim

k'

kk

k1k

k=α⇒

αα

−α=α⇒⎟⎟⎠

⎞⎜⎜⎝

⎛−=

∞→+

∞→

Deci rezultă că: ( )( ) ( ) ( ) 0ff;0xfecuatieiradacina =α=α=−α=α

Pentru a spori convergenţa iterativă a soluţiei, pentru a creşte precizia, ori pentru a reduce efortul de calcul din metoda iniţială a lui Newton au fost deduse alte forme îmbunătăţite sau simplificate.

Metoda Newton simplificată

( )( )0

'1 xfxf

xx nnn −=+

În cazul în care evaluarea derivatei funcţiei, în fiecare nou punct de aproximare, este costisitoare, formula lui Newton poate fi adaptată, în sensul reţinerii valorii calculate a derivatei în primul punct de aproximare pe parcursul unui număr de iteraţii, şi, eventual, schimbarea acestei valori a derivatei prin recalcularea într-un nou punct de aproximare cu :

Se reduce simţitor efortul computaţional, cu dezavantajul reducerii concomitente a convergenţei; sunt necesare mai multe iteraţii până la obţinerea unei soluţii dorite.

Metoda lui Halley

Din dezvoltarea în serie Taylor, dacă se reţine şi cel de-al doilea termen, se deduce o formulă de aproximare a soluţiei, care prezintă convergenţă ridicată, cu preţul unui efort de calcul de luat în seamă, datorită evaluării la fiecare iteraţie atât a primei derivate cât şi a celei de-a doua:

( )( ) ( )[ ] ( ) ( )nnnn

nnn

xfxfxfxf

xfxx

''2''1

2

2

⋅⋅−±

⋅−=+

Metoda secanteiDacă se urmăreşte eliminarea evaluării derivatei, aceasta se aproximează cu formula de mai jos, ştiut fiind că ecuaţia dreptei care trece prin două puncte ale funcţiei, este o dreaptă care poate înlocui tangenta la funcţie:

( ) ( ) ( )1

1'

−

−

−−

=nn

nnn xx

xfxfxf

care se introduce în expresia iniţială a lui Newton:

( ) ( ) ( )⎥⎦⎤

⎢⎣

⎡−−

⋅−=−

−+

1

11

nn

nnnnn xfxf

xxxfxx

Pentru aplicarea acestei metode trebuie cunoscute primele două aproximaţii iniţiale x0 şi x1!!!

Metoda lui Newton este o metodă locală, în sensul că procesul iterativ converge doar dacă aproximaţia iniţială este aleasă suficient de aproape de valoarea adevărată.

Convergenţa metodei lui Newton şi a variantelor alternative este mai rapidă decât a metodei bisecţiei. S-a observat că soluţia s-a identificat cu o precizie satisfăcătoare în primele trei iteraţii.

În alegerea uneia sau alteia dintre metodele numerice de soluţionare a ecuaţiilor trebuie să se estimeze efortul de calcul implicat pentru atingerea unei precizii, dificultatea de evaluare repetată a funcţiei şi a derivatelor ei.

Concluzii asupra metodei lui Newton

⎝ ⎠

f x( ) x2( ) x−sin x( )+ ex5

−

1

x

x2

x

x2

1

x2

1

x

⎛⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎠

+:= \\ functia atasata:

Reprezentarea grafica in vederea selectarii intervalului in care se afla solutia, sau solutiile:

x 0 0.01, 1..:=

0 0.2 0.4 0.6 0.8 1

2

2

f x( )

x

Metoda Newton-Raphson (metoda tangentei)

\\ ecuatia a carei solutii pozitive trebuie identificate;x2( ) x−

sin x( )+ ex5−

1

x

x2

x

x2

1

x2

1

x

⎛⎜⎜⎜⎜⎝

⎞⎟⎟⎟⎟⎠

+ 0

Algoritmi MathCad

--> functie Rolle f a( ) f b( )⋅ 0< f a( ) f b( )⋅ 0.037−= se verifica

--> derivata functiei f' x( )xf x( )d

d:= \\ trebuie sa fie pozitiva pe intervalul

considerat

Conditiile de unicitate a solutiei pe un interval selectat:

a 0.1:= b 0.2:= x a a 0.001+, b..:=

Se urmareste determinarea radacinii din intervalul [0;0.2], procedeul de identificare a celei din intervalul [0.8;1] fiind analog.

0.1 0.12 0.14 0.16 0.18 0.22

4

6

f' x( )

x

--> alegerea aproximatiei initiale: x0 0.15:= (din grafic)

f'' x( ) 2xf x( )d

d

2:= f x0( ) f'' x0( )⋅ 0≥ f x0( ) f'' x0( )⋅ 1.258= \\ se verifica;

--> stabilirea unui numar de iteratii care sa verifice toleranta admisa: N 4:= k 2 N..:=

--> impunerea unei erori de calcul: ε 10 3−:=

--> algoritmul de iterare:

g x( )f x( )f' x( )

:= xk xk 1−

f xk 1−( )f' xk 1−( )−:= xN 0.168630073370396= \\ rezultatul

numeric cu 15 zecimale

f xN( ) xN xN 1−−+ε

5≤ f xN( ) xN xN 1−−+ 0.000000000197321=

\\ relatia de verificare a nedepasirii erorii impuse.

ε

50.0002=

4. O pila de tensiune electromotoare imprimata U.e si de rezistenta interioara R.i alimenteaza un circuit electric format din rezistentele reglabile R.1 si R.2, conectate in tandem conform figurii de simbolizare. La deplasarea cursorului comun, rezistorul R.1 isi modifica rezistenta liniar cu distanta, iar rezistorul R.2 cosinusoidal cu miscarea rectilinie a cursorului. Sa se determine cate unitati de masura trebuie sa se deplaseze cursorul, din cursa totala, astfel incat sa se ajunga la o valoare nula a intensitatii curentului reglat prin rezistorul cu variatie liniara.

ΩR1 300=R1 R1 d( ):=cm d 10:=--> cursa totala

R2 x( ) B cos x( )⋅:=R1 x( ) A x⋅:=

Ω

cmB 25:=

Ω

cmA 30:=

AIr 0:=VUe 24:=ΩRi 0.15:=

I Ir I2+

Ri I⋅ R1 x( ) Ir⋅+ R1 R1 x( )−( ) I⋅+ Ue IrUe R2 x( )⋅

Ri R1+ R1 x( )−( ) R2 x( ) R1 x( )+( )⋅ R1 x( ) R2 x( )⋅+ <=>

R1 x( )− Ir⋅ R2 x( ) I2⋅+ 0

Modelul matematic adoptat pentru rezolvarea cerintei impuse porneste de la considerarea teoremelor lui Kirchhoff, care aplicate pe circuit, formeaza un sistem liniar de ecuatii, a carui solutie de interes se exprima analitic cu regula lui Cramer. Se figureaza etapele modelului analitic:

Dat fiindca s-a ajuns la o ecuatie neliniara avand ca necunoscuta deplasarea rectilinie x, se apeleaza la o metoda numerica de solutionare a ecuatiilor neliniare. Metoda destinata acestei probleme, datorita largului uz in electrotehnica, este Newton-Raphson, cu corespondent in geometrie: tangenta repetata la curba functiei care se ataseaza ecuatiei. Procesul iterativ Newton-Raphson de ordinul 2 presupune aproximarea initiala, succesiv, a doua solutii pentru ecuatie.

Domeniul de definitie al functiei se considera spatiul total posibil de parcurs al cursorului, iar pentru acesta se verifica eligibilitatea conditiilor de aplicare a metodei.

f x( ) IrA Ue⋅ cos x( )⋅

Ri R1+( ) A cos x( )⋅ B x⋅+( )⋅ B2 x2⋅−

−:= Domeniul de definitie: [0;d]

(afisarea numerica a coeficientilor functiei, cu ajutorul operatorului de evaluare simbolica definit in Mathcad)

f x( ) 720−cos x( )

9004.50cos x( )⋅ 7503.75x⋅ 625 x2⋅−+( )

⋅→

Se reprezinta grafic functia descrisa in vederea stabilirii domeniului restrans de definire:

0 1.25 2.5 3.75 5 6.25 7.5 8.75 10

0.2

0.2

f x( )

x

Se observa ca la valoarea nula a intensitatii curentului electric se ajunge in mai multe pozitii ale cursorului. Prima pozitie a cursorului in care curentul prin latura de circuit se anuleaza este situata, dupa reprezentarea grafica, in intervalul [0;3] mm.

Conditiile de unicitate ale solutiei pe intervalul selectat:

a 0:= b 3:= x a a 0.01+, b..:= (esantionarea intervalului pentru sporirea acuratetei reprezentarii)

--> functie Rolle f a( ) f b( )⋅ 0< 1= (atestarea booleana a verificarii conditiei)

--> derivata functiei f' x( )xf x( )d

d:= --> trebuie sa fie pozitiva pe intervalu

considerat

0 1 2 3

0.2f' x( )

x

--> alegerea aproximatiei initiale: x0 1.58:= (din grafic)

f'' x( ) 2xf x( )d

d

2:= f x0( ) f'' x0( )⋅ 0≥ 1= x1

a f x0( )⋅ x0 f a( )⋅−

f x0( ) f a( )−:= x1 1.567=

--> stabilirea unui numar de iteratii: N 10:= k 2 N..:=

0 1.25 2.5 3.75 5 6.25 7.5 8.75 10

0.2

0.2

f x( )

x

--> algoritmul de iterare:

xk xk 1−

f xk 1−( )f' xk 1−( )−:= xN 1.571= cm

x 1.58:= (aproximare initiala a solutiei)

Given

IrA Ue⋅ cos x( )⋅

Ri R1+( ) A cos x( )⋅ B x⋅+( )⋅ B2 x2⋅−

− 0

sol Find x( ):= sol 1.570796= cm

(punctul de pe ruta de deplasare a cursorului in care se anuleaza valoarea intensitatii curentului electric din portiunea de circuit cu rezistor cu variatie liniara).

Programul Mathcad ofera posibilitatea solutionarii rapide a ecuatiei deduse in aceasta problema, prin apelul blocului de functii Given- Find. Mai jos se exemplifica acest mod de rezolvare, urmand ca apoi sa se compare cele doua rezultate obtinute pentru cerinta problemei, in conditiile in care se va considera ca referinta afisajul numeric al blocului definit in Mathcad. Se mentioneaza, totusi, ca in spatele rulajului predefinit in programul Mathcad se afla o metoda numerica de tipul celei etapizate, cu precizie impusa.

Eroare sol rez−:= Eroare 8.57 10 9−×=

Desi problema se incadreaza in tendinta generala de abordare teoretica a posibilitatilor de aplicare a meastfel de situatie, precum cea propusa in problema isi poate gasi aplicare in circuite de pozitionare precivariatia curentului. Precizia de aproximare a pozitiei nu trebuie sa depaseasca 3 zecimale.