Newton – módszer A húrmódszernél és a szelőmódszernél az F(x) függvény gyökének közelítéséhez a függvény húrját használtuk. Hatásosabb a módszer akkor, ha érintőkkel dolgozunk. Def.: Legyen x0 az F(x) = 0 egyenlet x* gyökének elég jó közelítése. Húzzuk meg az F(x) függvény (x0, F(x0)) pontbeli érintőjét, és jelölje x1 az érintő x tengellyel vett metszéspontját. Ezután húzzuk meg az F(x) függvény (x1, F(x1)) pontbeli érintőjét és jelölje x2 az érintő x tengellyel vett metszéspontját… Íly módon egy x0, x1, x2,… sorozatot nyerünk. Xn+1 = xn – F(xn)/F’(xn)
(n = 0, 1, 2,…)
Tétel: Tegyük fel, hogy x0 és x* között F(x) kétszeresen differenciálható, F’(x)≠0, F’’(x)≠0, valamint F(x0)*F’’(x0) > 0. Ekkor a Newton módszer konvergens. Ha az x* gyököt a húr és a szelő-módszerekhez hasonlóan valamely intervallumba szorítjuk, akkor az intervallum azt a végpontját kell kezdeti közelítésnek választanunk, amelyre: F(x0)*F’’(x0)>0. A Taylor-formula alapján k>=0 esetén azt nyerjük, hogy: 0=F(x*)=F(xk)+F’(xk)*(x*-xk)+(1/2)*F’’(ξ)*(x*-xk)2 Valamint a: Xk+1 = xk – F(xk)/F’(xk) formulát átrendezve: 0=F(xk)+F’(xk)*(xk+1-xk) adódik, amelyet kivonva a Taylor-formulából nyert egyenlőségből azt kapjuk hogy: 0=F’(xk)*(x*-xk+1)+(1/2)*F’’(ξ)*(x*-xk)2 formulát kapjuk. Feltéve, hogy az x* értékét és az (xk) sorozatot tartalmazó valamilyen intervallumba 0<m<=|F’(x)| és M>=|F’’(x)|,akkor az előbbi formula: |ξk+1|<=M/2m*|ξk|2 alakban írható fel. Mindkét oldalt megszorozva K=M/2m-el a dk=K*ξk mennyiségekre a dk+1<=dk2 összefüggéshez jutunk. Feltéve, hogy x0 x* hoz olyan közel van, hogy d=d0<1, akkor teljes indukcióval könnyű belátni, hogy k>=0 esetén: dk+1<=d2k+1
Ez az egyenlőtlenség a Newton Módszer hibabecslő formulája.
Példaprogram: x*x*x-3*x*x-x+9 import math, os from abrazol import * rajzol(["(x*x*x)-3*(x*x)-x+9"],-5,5,60,20,["red"]) # a fuggveny def func1(x): return (x*x*x)-3*(x*x)-x+9; #elso derivalt def func2(x): return 3*x*x-6*x-1; #beolvasas x0=input("Adja meg az intervallum kezdetet: ")*1.0 eps=input("Adja meg a hibakorlatot: ")*1.0 nevezo=func2(x0) szamlalo=func1(x0) xe=x0-(szamlalo/nevezo) sz=1 while (abs(x0-xe)>eps): x0=xe szamlalo=func1(x0) nevezo=func2(x0) print sz print xe xe=x0-(szamlalo/nevezo) sz=sz+1 Hibabecslés: F(x)=x3-3x2-x+9 |F’(x)|=|3x2-6x-1|>=8=m |F’’(x)|=|6x-6|<=18=M K=M/2m K=18/16=9/8 ξk=xk-x* dk=K*ξk d=(K*ξ0)2k+1 x0=-2 x1=-1,609 ξ0=x0-x1=-0.391 d0=(9/8)*-0,0391=-0,0439875 dk+1<=d2k+1 d4<=(-0,0439875)9 9/8*|ξ5|<(-0,0439875)9 |x4-x*|<8/9*(-0,0439875)9
Futási eredmény: /usr/bin/python -u "/home/gajdosr/Python/oldal/Newton.py" Adja meg az intervallum kezdetet: -2 Adja meg a hibakorlatot: 0.000001 1 -1.60869565217 2 -1.52839805339 3 -1.52510768074 4 -1.52510225483
Módosított Newton-módszer: A Newton-módszer esetében újabb közelítés számításakor f és f’ egy-egy függvényértékét kell kiszámolnunk! Az intervallum felezésnél, a húr- és a szelőmódszernél lépésenként csak f(xk) értékét kell számolnunk! Így a Newton-módszer műveletigénye nagyobb, mint a többi módszeré! Ha f’ a gyök környezetében alig változik, akkor nem vétünk hibát, ha a Newton formulába f’(xk) helyett f’(x0) értékét írjuk. Tehát e formulát írjuk át: Erre:
Az így kapott formulát módosított Newton-módszernek nevezzük. Ezzel a módszer műveletigényét jelentősen lecsökkentjük, hiszen f’ értékét csak az iteráció megkezdésekor kell kiszámolnunk. Példa: Az x^3-3x^2-x+9 = 0 egyenlet valós gyökét a módosított Newton-módszerrel határozzuk meg. x0 = -2 választással f’ (x0) = 23 . 13 lépést kell számolnunk a gyök 6 tizedesre való meghatározásához. A számítás eredményeit a táblázat ismerteti. import math, os from abrazol import * rajzol(["(x*x*x)-3*(x*x)-x+9"],-5,5,60,20,["red"]) def func1(x): return (x*x*x)-3*(x*x)-x+9; def func2(x): return 3*x*x-6*x-1; x0=input("Kerem az intervallum kezdetet: ")*1.0 eps=input("Kerem a hibakorlatot: ")*1.0 nevezo=func2(x0) szamlalo=func1(x0) xe=x0-(szamlalo/nevezo) lepessz=1 while (abs(x0-xe)>eps): x0=xe szamlalo=func1(x0) print lepessz print xe xe=x0-(szamlalo/nevezo) lepessz=lepessz+1
Forrás: Szidarovszky Ferenc – Bevezetés a numerikus módszerekbe
Lagrange interpoláció: import math, os from abrazol import * n=input("Adja meg az alappontok szamat:") i=0 xtomb=[ None ] * n ytomb=[ None ] * n tomb=[ None ]*(n+1) while(i
Polinom interpoláció Függvényközelítések Azzal a kérdéssel foglalkozik, hogy a diszkrét pontokban adott függvényekhez hogyan lehet jól kezelhető, az adott pontokra minél jobban illeszkedő függvényeket konstruálni. A legkönnyebben kezelhető és a legkedvezőbb analitikus tulajdonságokat követő függvények a véges fokszámú polinomok, így a gyakorlati esetek nagy részében polinom közelítésekkel dolgozunk. Az adott pontokra való jó illeszkedésük szempontjából a polinomokkal való közelítések három típusát különböztetjük meg: 1. Interpoláció 2. A legkisebb négyzetek módszere 3. Csebisev-féle közelítés Az interpolációs polinomok az alappontokban ugyanazokat az értékeket veszik fel, mint az adott függvény, a legkisebb négyzetek és a Csebisev-féle közelítés módszerével nyert polinomok az alappontokban az adott függvényértékeknek csak közelítését adják. Interpoláció Az y=f(x) függvény értékkészlete legyen ismert az x0, x1, …, xn pontsorozaton, azaz y0=f(x0), y1=f(x1), …, yn=f(xn) Az x0, x1, …, xn pontsorozatot a továbbiakban interpolációs alappontoknak nevezzük. Az interpoláció célja, hogy olyan függvényt határozzunk meg, amely az [x0; xn] intervallumban közelítőleg megadja az alappontoktól eltérő helyeken is a függvényértékeket. Az eljárás lényege az, hogy az f(x) függvényt olyan F(x) függvénnyel közelítjük, amely az (xi;yi) (i=0,1, …, n) pontokban, az ún. kollokációs pontokban megegyezik f(x)-szel, azaz F(xi)=f(xi)≡yi (i=0,1, …, n). Az F(x) függvény előállítására szolgáló eljárást interpolációnak, az F(x) függvényt pedig interpolációs függvénynek nevezzük. Az F(x) függvény p(x)-szel jelölt polinom. Polinom interpoláció A polinom interpoláció a lineáris interpoláció egy általános fajtája. Egy y=p(x) polinom meghatározását jelenti, mely keresztül megy a (x1,y1), (x2,y2),…,(xn,yn) pontokon. Tehát adott n db pont ahol az egyes xi értékek mind különbözők, minden p(xi)=yi, (i eleme 0..n) és a polinom fokszáma legfeljebb n-1 lesz. A keresett p(x) polinomra minden esetben teljesülni kell a következő feltételeknek: 1. xi!=0, i eleme 0..n, 2. xi=xj akkor i=j, 3. a következő mátrix determinánsa nem 0 1, x0, x02,…, x0n-1 1, x1, x12,…, x1n-1 … 1, xn, xn2,…, xnn-1
Ha ezek teljesülnek a következő egyenletrendszert kell megoldani (ez az ún. Vandermonde mátrix):
A keresett polinom a következő alakban írható fel az ai-k ismeretében: import math, sys def mxprint(m): for i in range(size): for j in range(size): print m[i][j], print "" def mkmatrix(rows, cols): mk = [ None ] * rows for i in range(rows): mk[i] = [0] * cols for j in range(cols): mk[i][j] = 0 return mk def delete(mx,sor,oszlop): m=mkmatrix(len(mx)-1,len(mx)-1) sorindex=-1 for i in range(len(mx)): if(i!=sor): sorindex+=1 oszlopindex=-1 for j in range(len(mx)): if(j!=oszlop): oszlopindex+=1 m[sorindex][oszlopindex]=mx[i][j] return m def mxdet(m): ejel=-1 ret=0 for i in range(len(m)): if (len(m)==2): ret=(m[0][0]*m[1][1])-(m[1][0]*m[0][1]) else: ret=ret+(-1)*ejel*m[0][i]*(mxdet(delete(m,0,i))); ejel=-ejel return ret def mod(m1,m2,el): mke=mkmatrix(len(m2),len(m2)) for i in range(len(mke)): for j in range(len(mke)): mke[i][j]=m1[i][j] mke[i][el]=m2[i] return mke def cramer(m1,m2): xi=mkmatrix(len(m2),len(m2)) for i in range(len(m2)): xi[i]=mxdet(mod(m1,m2,i))/(mxdet(m1)*1.0) return xi
def fuggveny(x): i=1 ertek=p[0] while(i
Az intervallumfelezés módszere import math def f(a): return a*a*a*a*a+2*a*a-0.5 def intfel(a,b): c=(a+b)/2.0 if ((abs(a-b)>1e-11)and(f(c)!=0)): if (f(c)*f(b)<0): return intfel(c,b) elif (f(c)*f(a)<0): return intfel(a,c) else: return c def intfel1(a,b): d=0 c=(a+b)/2.0 while((abs(a-b)>1e-11) and (f(c)!=0)): c=(a+b)/2.0 if (f(c)*f(b)<0): a=c elif (f(c)*f(a)<0): b=c d+=1 return d a=0 b=2 if (f(a)*f(b)<0): print "a kozelito megoldas :",intfel(a,b) print intfel1(a,b),"lepesben oldotta meg" else: print "Nemjo a megadott intervallum"
A legkisebb négyzetek módszere Eddig már két függvényközelítési módszerrel foglalkoztunk, a Lagrange polinomokkal és a Taylor polinomokkal. A Lagrange polinomoknál minden alappontban egy mérési eredményünk van (ami lehet valódi mérés eredménye, de lehet kiszámított függvényérték is) és megköveteljük azt, hogy a függvényt közelít? polinom a megadott alappontban a megadott értéket vegye fel. A Taylor polinomok esetében egy pontban a deriváltak értékét adjuk meg (illetve mérjük, ha ilyen mérést meg tudunk valósítani) és olyan polinomot konstruálunk, amelynek deriváltjai az adott pontban a megadott derivált értékek. A legkisebb négyzetek módszere a fenti módszerek egy általánosítása, ugyanis a gyakorlatban meg kell engednünk azt is, hogy egy függvényérték meghatározására több mérést is végezhessünk. Ekkor azonban nem köthetjük ki, hogy a közelítő függvény milyen értéket vegyen fel, hiszen a mérési eredmények rendszerint nem azonosak, így nincs is megadott függvényérték. A másik általánosítás abban lehetséges, hogy nem kell ragaszkodni a polinomokhoz, szinte minden függvényfajta előfordulhat illesztő függvényként.
A legkisebb négyzetek módszerének általános megfogalmazása Tegyük fel, hogy egy f ( x, a1 , a 2 ,..., a m ) egyelőre ismeretlen függvény értékére az x1 , x 2 ,...x s alappontokban méréseket végzünk. Ennek eredményeként az y j = f ( x j , a1 , a 2 ,...a n ) j = 1, 2, ..., s értékekre kapjuk az
, ...,
,
, ...,
,
........................................
, ...,
,
mérési eredményeket. ( s1 , s 2 ,.., s k ) nem feltétlenül egyenlőek, vagyis nem minden pontban kell ugyanannyi mérést végezni (de lehet). A fő feladat az f ( x, a1 , a 2 ,..., a m ) függvényben, amelynek alakja 2 ( m − 1) adott (pl. egy polinom, f ( x, a1 , a 2 ,..., a m ) = a1 + a 2 x + a3 x + ... + a m x ) az a1 , a 2 ,...a m határozatlan együtthatók értékének meghatározása úgy, hogy az így kapott f függvény értékének eltérése a mérési értékektől az alappontokban a lehető legkisebb legyen. Az eltérést a függvényérték és a mérési értékek különbségének négyzetével mérjük. Így a kapott feladat egy többváltozós függvény szélsőértékének meghatározása. A négyzetes eltérést megadó függvény a következő: f ( x, a1 , a 2 ,..., a m ) =
.
Keresendő tehát az f ( x, a1 , a 2 ,..., a m ) függvény minimumhelye, ahol változók az a1 , a 2 ,...a m paraméterek. A többváltozós függvények elméletéből tudjuk, hogy ott lehetnek szélsőérték helyek, ahol a függvény első parciális deriváltjai eltűnnek. Esetünkben ez a következő egyenletek teljesülését jelenti:
=
= ...
= A szélsőérték létezésének elégséges feltételeivel ilyen általánosan nem foglalkozunk. Abban az esetben, ha a meghatározandó függvény polinom vagy olyan függvény, amelyben az ismeretlen a j , j = 1 , ..., m paraméterek lineárisan fordulnak elő, lineáris legkisebb négyzetek módszeréről beszélünk. Ennek speciális esete, amellyel külön is foglalkozunk az, amikor f alakja f ( x, a1 , a 2 ) = a1 + a 2 x vagyis a lineáris függvény, amelyet a statisztikában lineáris regressziónak neveznek. De ugyanebben az értelemben beszélhetünk parabolikus, harmadfokú, ... stb. regresszióról is, ha az illesztésre használt függvény parabola, harmadfokú polinom, stb. Minden ilyen esetben a fenti megoldandó egyenletrendszer lineáris egyenletrendszer lesz. Nemlineáris regresszióról akkor beszélünk, ha az illesztendő függvény a meghatározandó paramétereket nemlineárisan tartalmazza. Ekkor a szükséges feltételeket megfogalmazó egyenletrendszer nemlineáris egyenletrendszer lesz. ### lnm.py # -*- coding: iso-8859-1 -*import math, os, sys n=input("az alappontok szßma:") i=0 xtomb=[ None ] * n ytomb=[ None ] * n tomb=[ None ]* n while(i
while(i
k=0 q=[ [ None ] * n] * n p=[ [ None ] * n] * n while(k
while(l
Deteminans: from random import randrange def mxprint(m): for i in range(len(m)): for j in range(len(m)): print m[i][j], print "" def mkmatrix(rows, cols): mk = [ None ] * rows for i in range(rows): mk[i] = [0] * cols for j in range(cols): mk[i][j] = 0 return mk def mkrandommatrix(rows,cols): mk = [ None ] * rows for i in range(rows): mk[i] = [0] * cols for j in range(cols): mk[i][j] = randrange(20) return mk def mxdet(m): ejel=-1 ret=0 for i in range(len(m)): if (len(m)==2): ret=(m[0][0]*m[1][1])-(m[1][0]*m[0][1]) else: ret=ret+(-1)*ejel*m[0][i]*(mxdet(delete(m,0,i))); ejel=-ejel return ret def delete(mx,sor,oszlop): m=mkmatrix(len(mx)-1,len(mx)-1) sorindex=-1 for i in range(len(mx)): if(i!=sor): sorindex+=1 oszlopindex=-1 for j in range(len(mx)): if(j!=oszlop): oszlopindex+=1 m[sorindex][oszlopindex]=mx[i][j] return m def main(): #matrix=([1,3,2,2],[-2,6,2,6],[3,6,2,5],[1,2,1,1]) #matrix=([-2,1,2,4,-1],[3,1,1,-4,5],[-6,6,7,6,11],[11,10,-13,-9,6],[3,5,5,3,-8]) matrix=([-2,1,2,4,-1,5],[3,1,1,-4,5,-6],[-6,6,7,6,11,7],[11,10,-13,-9,6,8],[3,-5,5,3,-8,9],[2,2,11,3,-4,3]) #matrix=mkrandommatrix(11,11) # print matrix print mxdet(matrix) main() #by:TGT
Inverz def mxprint(m): for i in range(len(m)): print "" for j in range(len(m)): print m[i][j] def mkmatrix(rows, cols): mk = [ None ] * rows for i in range(rows): mk[i] = [0] * cols for j in range(cols): mk[i][j] = 0 return mk def mxdet(m): ejel=-1 ret=0 for i in range(len(m)): if (len(m)==2): ret=(m[0][0]*m[1][1])-(m[1][0]*m[0][1]) else: ret=ret+(-1)*ejel*m[0][i]*(mxdet(delete(m,0,i))); ejel=-ejel return ret def delete(mx,sor,oszlop): m=mkmatrix(len(mx)-1,len(mx)-1) sorindex=-1 for i in range(len(mx)): if(i!=sor): sorindex+=1 oszlopindex=-1 for j in range(len(mx)): if(j!=oszlop): oszlopindex+=1 m[sorindex][oszlopindex]=mx[i][j] return m def invertal(m): m1=mkmatrix(len(m),len(m)) ejel=-1 for i in range(len(m)): for j in range(len(m)): m1[j][i]=(-1)*ejel*mxdet(delete(m,i,j)) ejel=-ejel for i in range(len(m1)): for j in range(len(m1)): m1[i][j]=m1[i][j]*1.00 m1[i][j]=m1[i][j]/mxdet(m) return m1 def main(): #matrix=([1,-1,2],[2,-1,3],[1,-2,4]) #matrix=([1,2,3],[1,4,0],[-1,1,-1]) matrix=([-2,1,2,4,-1],[3,1,1,-4,5],[-6,6,7,6,11],[11,10,-13,-9,6],[3,-5,5,3,8]) mxprint(invertal(matrix)) main() #by:TGT
Lineáris egyenletrendszer: def mxprint(m): for i in range(len(m)): print "" for j in range(len(m)): print m[i][j] def mkmatrix(rows, cols): mk = [ None ] * rows for i in range(rows): mk[i] = [0] * cols for j in range(cols): mk[i][j] = 0 return mk def mxdet(m): ejel=-1 ret=0 for i in range(len(m)): if (len(m)==2): ret=(m[0][0]*m[1][1])-(m[1][0]*m[0][1]) else: ret=ret+(-1)*ejel*m[0][i]*(mxdet(delete(m,0,i))); ejel=-ejel return ret def delete(mx,sor,oszlop): m=mkmatrix(len(mx)-1,len(mx)-1) sorindex=-1 for i in range(len(mx)): if(i!=sor): sorindex+=1 oszlopindex=-1 for j in range(len(mx)): if(j!=oszlop): oszlopindex+=1 m[sorindex][oszlopindex]=mx[i][j] return m def mod(m1,m2,el): mke=mkmatrix(len(m2),len(m2)) for i in range(len(mke)): for j in range(len(mke)): mke[i][j]=m1[i][j] mke[i][el]=m2[i] return mke def cramer(m1,m2): xi=mkmatrix(len(m2),len(m2)) for i in range(len(m2)): xi[i]=mxdet(mod(m1,m2,i))/(mxdet(m1)*1.0) return xi def main(): matrix1=([3,2,1],[5,0,3],[9,4,3]) matrix2=([1,2,3]) #matrix1=([1,3,2,2],[-2,6,2,6],[3,6,2,5],[1,2,1,1]) #matrix2=([1,2,3,4]) print cramer(matrix1,matrix2) main()
