Seminar Numericne metode Mathematica

Univerza v Ljubljani FAKULTETA ZA STROJNIŠTVO NUMERIČNE METODE 2011/12 Seminar: Numerične metode PD-MK Datoteka: Seminar_numericne_metode_5, Zadnja sprememba: 17.09.13 Žiga Gosar 17.09.13 Seminar: Numerične metode KAZALO 1 Naloga 1 ............................................................................................................................................................ 4 1.1 Kompozitno Simpsonovo pravilo ................................................................................................................... 4 1.2 Definicija naloge (str. 246, naloga 5) ............................................................................................................. 4 1.3 Programska koda – Mathematica ................................................................................................................... 4 1.4 Rešitev fizikalnega problema ......................................................................................................................... 7 2 Naloga 2 ............................................................................................................................................................ 9 2.1 Taylorjeva metoda višjega reda ...................................................................................................................... 9 2.2 Definicija naloge (str. 272, naloga 6d) ......................................................................................................... 10 2.3 Programska koda – Mathematica ................................................................................................................. 10 2.3.1 Kubična Hermitova interpolacija ................................................................................................................. 12 3 Naloga 3 .......................................................................................................................................................... 14 3.1 Adaptivna ocena koraka, kontrola napake in Runge-Kutta-Fehlberg-ova metoda .......................................... 14 3.2 Definicija naloge (str. 288, naloga 3d) ......................................................................................................... 15 3.3 Programska koda – Mathematica ................................................................................................................. 15 3.3.1.1 4 Rešitev........................................................................................................................................................... 17 Naloga 4 .......................................................................................................................................................... 18 4.1 Numerično reševanje diferencialnih enačb ................................................................................................... 18 4.2 Definicija naloge (str 323, naloga 6) ............................................................................................................ 19 4.3 Programska koda – Mathematica ................................................................................................................. 19 4.3.1.1 5 Rešitev........................................................................................................................................................... 20 Naloga 5 .......................................................................................................................................................... 23 5.1 Teorija aproksimacij in Čebišev polinom...................................................................................................... 23 5.2 Definicija naloge (str. 516, naloga 1c) .......................................................................................................... 24 5.3 Programska koda – Mathematica ................................................................................................................. 24 5.3.1.1 6 Rešitev........................................................................................................................................................... 25 Naloga 6 .......................................................................................................................................................... 27 6.1 Teorija paproksimacij in Fast Foutier Transformacija................................................................................... 27 6.2 Definicija naloge (str. 547, naloga 5c) .......................................................................................................... 28 6.3 Programska koda – Mathematica ................................................................................................................. 28 7 6.3.1.1 Rešitev:.......................................................................................................................................................... 30 6.3.1.2 Rešitev........................................................................................................................................................... 31 Naloga 7 .......................................................................................................................................................... 33 7.1 QR algoritem ............................................................................................................................................... 33 7.1.1 QR iteracija ................................................................................................................................................. 33 7.2 Definicija naloge (str. 595, naloga 2d) ......................................................................................................... 34 7.3 Programska koda – Mathematica ................................................................................................................. 34 7.3.1.1 8 Rešitev........................................................................................................................................................... 35 Naloga 8 .......................................................................................................................................................... 36 Žiga Gosar Stran 2/49 Seminar: Numerične metode 8.1 Kontinuitetna metoda ................................................................................................................................... 36 8.2 Definicija naloge (str. 642, naloga 2) ........................................................................................................... 38 8.3 Programska koda – Mathematica ................................................................................................................. 38 8.3.1.1 9 Rešitev........................................................................................................................................................... 40 Naloga 9 .......................................................................................................................................................... 41 9.1 Hiperbolične parcialne diferencialne enačbe ................................................................................................. 41 9.2 Definicija naloge (str 725, naloga 7) ............................................................................................................ 42 9.2.1 Programska koda – Mathematica ................................................................................................................. 43 10 Primerjava numeričnih metod ....................................................................................................................... 47 11 Literatura ....................................................................................................................................................... 49 Žiga Gosar Stran 3/49 Seminar: Numerične metode 1 Naloga 1 1.1 Kompozitno Simpsonovo pravilo Simpsonovo pravilo dovolj dobro aproksimira integracijo krivulje na dovolj majhnem intervalu (intervalu, na katerem je krivulja relativno gladka). Tako lahko integral na takem intervalu aproksimiativno zapišemo kot:  ba   a b   f  xdx  6  f  a   4 f  2   f b  b a Na intervalu, na katerem se krivulja obnaša preveč oscilativno ali na nekaterih mestih ni diferenciabilna, Simpsonova aproksimacija utegne dati slab rezultat. Zato je potrebno Simpsonovo pravilo razširiti v kompozitno Simpsonovo pravilo, ki interval integracije razdeli na več podintervalov, na katerih se vrši aproksimativna integracija s Simpsonovim pravilom. Prispevke posameznih intervalov seštejemo, pa imamo ob sodem naravnem številu n, h= (b-a)/n in xj= a+ jh:  b a f  xdx  n 2 1 n 2 1  b  a 4  4 h 2 4   f a f x f  x2 j 1   f  b    h f        2j  3 180 j 1 j 1  kjer zadnji člen predstavlja napako aproksimacije. Kadar aproksimiramo neprave integrale z neskončnimi mejami osnovne oblike:   a 1 dx xp jih je potrebno najprej transformirati v integrale s končnimi mejami z vpeljavo substitucije t= 1/x. Tako dobimo:   a 1 1 1 a 2   dx t   dt  p  0 x t  kjer je izpeljana oblika primerna za uporabo kompozitnega Simpsonovega pravilo. Kjer vpeljava substitucije ni možna, je potrebno integral predhodno deliti na več integralov. 1.2 Definicija naloge (str. 246, naloga 5) Predpostavimo, da telo mase m potuje navpično navzgor, začenjši s površja zemlje. Če zanemarimo vse upore razen gravitacije je ubežna hitrost v podana kot v2  2 gR z2 dz , kjer je z   1 x R kjer je R = 6378.137 km radij zemlje in g = 9.81 m/s2 gravitacijski pospešek na zemljini površini. Predpostavi izhodno hitrost v. 1.3 Programska koda – Mathematica Simpson[a_,b_]:=Module[{},Return[(b-a)/6 (f[a]+4 f[(a+b)/2]+f[b])]]; Žiga Gosar Stran 4/49 Seminar: Numerične metode Adapt[a_,b_,tol_]:=Module[{c},c=(a+b)/2; Sab=Simpson[a,b]; Sac=Simpson[a,c]; Scb=Simpson[c,b]; If[Abs[Sab-Sac-Scb]<tol, Return[Sac+Scb], Return[Adapt[a,c,tol/2]+Adapt[c,b,tol/2]]];]; Adapt2[a_,b_,tol_]:=Module[{c},c=(a+b)/2; X=Append[X,c]; Sab=Simpson[a,b]; Sac=Simpson[a,c]; Scb=Simpson[c,b]; If[Abs[Sab-Sac-Scb]<tol, Print[ Print[ " "ab, f x , " x " " "ca, f "cb, f ", Sab]; x ," x ", Sac]; x ," x ", Scb ]; Print[ Print["Je v tolerančnem območju = ",tol]; Print["-----------"]; Return[Sac+Scb], Return[Adapt2[a,c,tol/2]+Adapt2[c,b,tol/2]]];]; f[x_]=1/x2; a=1.0; b=100.0; X={a,b}; val=Adapt2[a,b,10-6]; X=Sort[X]; Print["f[x] = ",f[x]]; Print[""]; " "b, f x , " x ", val a Print[ ]; Print[""]; Print["Točke na abcisi {xk} za interval [xk-1,xk] s Simpsonovim pravilom"]; xk Print[" x k 1 f[x]x","  ","(xk-xk-1)/6(f[xk-1]+4f[(xk-1+xk )/2]+f[xk])"]; Print["kjer so točke"]; Print["{xk} = ",X]; Print["\n"]; tol1=0.001; val1=Adapt[a,b,tol1] tol2=0.00001; val2=Adapt[a,b,tol2] tol3=0.0000001; val3=Adapt[a,b,tol3] Print[""] b 1 2 x val4= a x Print[""] Print["\n"]; Needs["Graphics`FilledPlot`"]; Needs["Graphics`Colors`"]; Žiga Gosar Stran 5/49 Seminar: Numerične metode X=Sort[X]; Y=f[X]; pts=Transpose[{X,Y}]; dots=ListPlot[pts,PlotStyle{Magenta, PointSize[0.01]}, DisplayFunctionIdentity]; gr=FilledPlot[f[x],{x,1,b},PlotRange{{1,10},{0,1}},PlotStyleMagenta, FillsCyan,AxesLabel{"x","y"}, DisplayFunction$DisplayFunction]; Show[gr,dots,DisplayFunction$DisplayFunction] Print["f[x] = ", f[x]]; Print["{xk} = ", X]; Print["Numerična vrednost integrala je ", val]; Rešitev: y 1.0 0.8 0.6 0.4 0.2 x 2 4 Numerična vrednost integrala je 6 8 10 0.99 Integrate[x^-2,{x,1,}] 1 Rezultat: v0  2 gR  0,990000002  11,06209474 km/s v0  2 gR  0,990000002  6,910171633 mi/s Rešitev v literaturi [1], str 775. 6.9450 mi/s Žiga Gosar Stran 6/49 Seminar: Numerične metode Rešitev fizikalnega problema 1.4 Na telo mase m, ki je na razdalji r (r  R – polmer Zemlje) od središča Zemlje deluje gravitacijska sila F   mM , r2 (χ je gravitacijska konstanta, M pa predstavlja maso Zemlje), ki da telesu pospešek enak: a  r   M . r2 To je diferencialna enačba gibanja telesa. r v r dr , dt dv dv dr dv   v. dt dr dt dr Ko vstavimo to v gornjo enačbo, sledi: v dv M   2 , dr r  vdv    M  r dr v2  M  C. 2 r 2 , Označimo z v0 začetno hitrost (hitrost izstrelitve), to je hitrost pri r = R in vstavimo v zadnjo enačbo; dobimo konstanto C, C v02  M  . R 2 Sledi: v2  M  v02  M     . r R 2 2   Ker telo zapusti zemljsko gravitacijsko polje, ima pri poljubno velikem r v>0. To pomeni, da mora biti izpolnjen pogoj: v02  M  . R 2 Od tu takoj dobimo iskano začetno hitrost 2 M . R v0  Ker je tečni pospešek g M 2 R  9,81 m/s2 , R  6,3 106 m, Žiga Gosar Stran 7/49 Seminar: Numerične metode sledi v0  2 gR  11,1178 v0  11, 2 km/s2 . Žiga Gosar Stran 8/49 Seminar: Numerične metode 2 Naloga 2 2.1 Taylorjeva metoda višjega reda Bistvo numeričnih metod je doseganje dovolj natančni aproksimacij k eksaknim rešitvam diferencialnih enačb z minimalnim trudom. Učinkovitost aproksimacije lahko ugoravljamo na več načinov. Eden izmed njih je sledenje lokalne napake, ki nam pove, kolikšno je odstopanje od rešitve diferencialne enačbe, pridobljene z diferenčno enačbo po vsakem koraku aproksimacije. Diferenčno enačbo zapišemo kot: wi 1  wi  h  ti , wi  , i  1,..., N 1 2.1 Njena lokalna napaka je:  i 1  h   yi 1  yi    ti , yi  , h i  1,..., N 1 Lokalna napaka Eulerjeve metode (ki je posebna oblika Taylorjeve metode z n = 1) pri i-tem koraku za začetni problem: y  f  t , y , je  i 1  h   a  t  b, yi 1  yi  f  ti , yi  , h ya    , 2.2 i  1,..., N 1 Kjer predstavlja yi = y(ti) eksaktno vrednost pri času ti. Napaka je lokalna, ker pomeni odstopanje pri določenem koraku ob predpostavki, da je rešitev ob prejšnjem koraku eksaktna. Kot taka je odvisna od diferencialne enačbe, velikosti koraka in od posameznega koraka aproksimacije. Ideja je torej izboljšati konvergenco rešitve z višanjem stopnje n Taylorjevega polinoma. Za začetni problem (enačba 2.2), ob predpostavki, da je rešitev diferencialne enačbe (n+1)-krat odvedljiva, lahko zapišemo Taylorjev polinom n-te stopnje za rešitev y(t) okoli ti za oceno pri ti+ 1 kot: y  ti 1   y  ti   hy  ti   h2 y  ti   2!  hn  n hn 1  n1 y  ti   y  i  n!  n  1! 2.3 y  t   f  t , y  t   Za nekatere ξi v intervalu (ti,ti+ 1). Zapišemo odvode: y  t   f   t , y  t   y k   t   f  k 1  t , y  t   Zapisane odvode zamenjamo v (enačba), pa imamo: y  ti 1   y  ti   hf  ti , y  ti    Žiga Gosar h2 f   ti , y  t i    2!  hn  n1 hn1 f ti , y  t i    f n  i , y  i    1 ! n! n    Stran 9/49 Seminar: Numerične metode Taylorjeva metoda stopnje n je diferenčna metoda, uporabljena na enačbi (2.3), pri čemer zanemarimo zadnji člen. Aproksimativno rešitev za (i+1)-vi korak torej zapišemo kot: wi 1  wi  hT  n  ti , wi  kjer je T 2.2 n  ti , wi   f  ti , wi   h f   ti , wi   2!  hn1  n1 f ti , wi  n! Definicija naloge (str. 272, naloga 6d) Podan je problem z začetnim pogojem y  1 y   y2 , 1  t  2 , t2 t 1 Z eksaktno rešitvijo y  t    : t y 1  1 c) Uporavi Taylorjevo metodo četrtega reda z h = 0.05 da oproksimirate rešitev in primerjajte rezultat z dejansko vrednostjo y. d) Uporabi rešitve generirane v delu naloge c) in uporabi kubično Hermitovo interpolacijo za aproksimiranje sledečih vrednosti y in jih primerjaj z dejanskimi vrednostmi. i. y 1.052  2.3 ii. y 1.555  iii. y 1.978  Programska koda – Mathematica (*Taylor 4.reda*) Clear[t,y,h,f]; Dy[t_,y_]=1/t^2-y/t-y^2 D1y[t_,y_]=D[Dy[t,y],t]+D[Dy[t,y],y]*Dy[t,y] (*T2 je do sem, brez D2y in D3y*) D2y[t_,y_]=D[D1y[t,y],t]+D[D1y[t,y],y]*Dy[t,y] D3y[t_,y_]=D[D2y[t,y],t]+D[D2y[t,y],y]*Dy[t,y] f[t_,y_,h_]=y+h (Dy[t,y]+D1y[t,y]h/2+D2y[t,y]h^2/6+D3y[t,y]h^3/24) tn=1; (*spremeni*) yn=-1; (*spremeni*) hn=0.05;(*spremeni*) Y={{tn,yn}}; For[i=0,i19,i++,AppendTo[Y,{tn+hn,f[tn,yn,hn]}]; {tn,yn}=Y[[-1]]] g1=ListPlot[Y]; g2=Plot[-1/t,{t,1,2}]; Show[g1,g2]; (*Računanje vrednosti*) t=1.0;(*spremeni za korak h*) h=0.05; y=w=-1;(*kopiraj rezultat*) y+h (1/t2-y/t-y2+1/2 h (-(2/t3)+y/t2+(-(1/t)-2 y) (1/t2-y/t-y2))+1/6 h2 (6/t4-(2 y)/t3+(-(1/t)-2 y) (-(2/t3)+y/t2)+(1/t2-y/t-y2)/t2+(1/t2-y/t-y2) (1/t2+(-(1/t)-2 y)2-2 (1/t2-y/t-y2)))+1/24 h3 (-(24/t5)+(6 y)/t4+(-(1/t)-2 y) (6/t4-(2 y)/t3)+(2 ((2/t3)+y/t2))/t2-(2 (1/t2-y/t-y2))/t3+(1/t2-y/t-y2) (-(2/t3)+(2 (-(1/t)-2 y))/t2-2 Žiga Gosar Stran 10/49 Seminar: Numerične metode (-(2/t3)+y/t2))+(-(2/t3)+y/t2) (1/t2+(-(1/t)-2 y)2-2 (1/t2-y/t-y2))+(1/t2-y/t-y2) (-(2/t3)+(2 (-(1/t)-2 y))/t2-2 (-(2/t3)+y/t2)-6 (-(1/t)-2 y) (1/t2-y/t-y2)+((1/t)-2 y) (1/t2+(-(1/t)-2 y)2-2 (1/t2-y/t-y2))))) Tabela 1: Taylorjeva metoda 4. reda t 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 T(4) -1.00000000 -0.95238125 -0.90909144 -0.86956594 -0.83333422 -0.80000103 -0.76923192 -0.74074199 -0.71428706 -0.68965661 -0.66666819 -0.64516289 -0.62500167 -0.60606235 -0.58823710 -0.57143044 -0.55555749 -0.54054254 -0.52631785 -0.51282263 -0.50000218 y(t) -1.00000000 -0.95238095 -0.90909091 -0.86956522 -0.83333333 -0.80000000 -0.76923077 -0.74074074 -0.71428571 -0.68965517 -0.66666667 -0.64516129 -0.62500000 -0.60606061 -0.58823529 -0.57142857 -0.55555556 -0.54054054 -0.52631579 -0.51282051 -0.50000000 napaka T(4) 0.00000000 0.00000030 0.00000053 0.00000073 0.00000089 0.00000103 0.00000115 0.00000125 0.00000135 0.00000144 0.00000152 0.00000160 0.00000167 0.00000174 0.00000181 0.00000187 0.00000194 0.00000200 0.00000206 0.00000212 0.00000218 Taylorjeva metoda 2. in 4. reda -0.50 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 -0.55 -0.60 -0.65 y(t) -0.70 T(4) T(2) -1/t -0.75 -0.80 -0.85 -0.90 -0.95 -1.00 t Diagram 1: Primerjava rezultatov Taylorjeve metode z eksaktno rešitvijo Žiga Gosar Stran 11/49 Seminar: Numerične metode 2.3.1 Kubična Hermitova interpolacija Naloga od nas zahteva, da izračunamo aproksimativne vrednosti v točkah, ki so pozicionirane med izračunanimi, y 1.052  , y 1.555 , y 1.978 . V primeru, da moramo določiti aproksimacijsko vrednost točke, ki je med točkama, lahko uporabimo Taylorjevo metodo četrtega ali drugega reda. Za določitev vrednosti med točkama uporabimo, na Taylorjevi metodi, linearno interpolacijo. Računamo aproksimacijo za točko y(1.052), ki je med točkama t = 1.05 in t = 1.1.  1.052  1.1   1.052  1.05  y 1.052      -0.95238125     -0.909144  0.95065  1.05  1.1   1.1  1.05  Eksaktna vrednost v točki y(1.052) = -0.950570, aproksimacija ima tako napako 0.000079. Tabela 2: Aproksimatvne vrednosti z uporabo linearne interpolcije Taylorjevega polinoma t 1,052 1,555 1,978 T(4) -0.950650 -0.643147 -0.505643 y(t) -0.950570 -0.643087 -0.505561 napaka T(4) 0.000079 0.000060 0.000082 Za boljšo aproksimacijo vrednosti v točki y(1.052), uporabimo kubično Hermitovo interpolacijo. To od nas zahteva aproksimacijo v točkah y'(1.05) ,y'(1.1) in aproksimacijo vrednosti v y(1.05) in y(1.1). y 1.05  in y 1.1  y 1.05 1 2   y 1.05  0.90702920 2 1.05 1.05 y 1.1 1 2   y 1.1  0.82644580 . 2 1.1 1.1 Podčtrane rezultate smo dobili s pomočjo podanih podatkov, ostale rezultate pa dobimo s pomočjo formule (glej tabelo ). Tabela 3: Računanje po kubični Hermitovi interpolaciji z z0  x0 f[z] Prve razlike   f  z0   f x0 Druge razlike Tretje razlike   f  z0 , z1   f  x0 z1  x0   f z ,z  f z ,z  f  z0 , z1 , z2    1 2   0 1  z2  z0 f  z1   f x0 f  z1 , z2   Žiga Gosar f  z2   f  z1  z2  z1   f  z0 , z1 , z2 , z3   f  z1 , z2 , z3   f  z0 , z1 , z2  z3  z0 Stran 12/49 Seminar: Numerične metode f  z2  f  x1  z2  x1       f  z1 , z2 , z3   f  z2 , z3   f  z1 , z2  z3  z1   f  z2 , z3   f  x1 f  z3   f x1 z3  x1 Tabela 4: Primer računanja kubične Hermitove interpolacije z 1.05 Hamilton T(4) f[z] -0.95238125 Prve odvod Druge odvod Tretje odvod Rezultat 0.75309168 -0.95057063 0.90702920 1.05 -0.95238125 -0.82466129 0.86579613 1.10 -0.90909144 1.10 -0.90909144 -0.78700671 0.82644580 y 1.052T (4)   0.95238   t  1.05 0.90702    t 1.05   0.82466    t 1.05   t 1.1 0.75309 Kubični Hermitov polinom ima tako obliko: y 1.052T (4)   0.95057063 2 2 Tabela 5: Aproksimativne vrednosti kubičnega Hermitovega polinoma t 1,052 1,555 1,978 Žiga Gosar Hermite(T(4)) -0.9505706 -0.6430884 -0.5055633 y(t) -0.9505703 -0.6430868 -0.5055612 napaka H(T(4)) 0.0000003 0.0000016 0.0000021 Stran 13/49 Seminar: Numerične metode 3 Naloga 3 3.1 Adaptivna ocena koraka, kontrola napake in Runge-Kutta-Fehlberg-ova metoda Aproksimacija je lahko bolj učinkovita, če uporabljamo spremenljivo velikost koraka, vendar je vpeljava take izboljšave lahko komplicirana. Vendar pa take aproksimacije ob spreminjanju velikosti koraka vključujejo tudi oceno lokalnih napak brez potrebe računanja njihovih aproksimacij v obliki višjih odvodov funkcije. Take aproksimacije niso uporabne le za aproksimativno integriranje, temveč tudi za aproksimativno reševanje začetnih problemov, saj omogočajo sprotno kontrolo napake. Diferenčna metoda: wi 1  wi  hi  ti , wi , hi  , i  0,1, , N 1 Ki aproksimira rešitev začetnega problema: y  f  t , y , a  t  b, ya    Pri neki prednastavljeni toleranci ε > 0 in minimalnem številu delitvenih točk zagotavlja, da globalna napaka |y(ti) - wi| nikoli ne preseže ε. Minimalno število delitvenih točk je pogojeno s spreminjanjem razdalje med njimi. Zapišemo aproksimacijo po Taylorjevi metodi reda n: wi 1  wi  hi  ti , wi , h  z lokalno napako  i 1  h   O  hn  in aproksimacijo po Taylorjevi metodi reda (n+1): wi 1  wi  hi  ti , wi , h  z lokalno napako  i 1  h   O  hn1  . Predpostavimo wi  y  ti   wi , lahko izpeljemo oceno lokalne napake aproksimacije reda n:  i 1  h   1  wi 1  wi 1  h Cilj je držati to lokalno napako v določenih mejah. To storimo z vpeljavo parametra K, ki je noedvisen od koraka h, tako da je:  i 1  h   Khn Pri čemer lahko pri novem koraku velikosti qh ob uporabi še obeh prvotnih aproksimacij n-tega in (n+1)-ega reda zapištemo:  i 1  qh   qn  wi 1  wi 1  h Sedaj lahko zapišemo odvisnost q in prednastavljene tolerance ε.  h q    wi 1  wi 1    1n Metoda Runge-Kutta-Fehlgerg uporablja to neenakost, pri čemer je za oceno lokalne napake 4. reda: Žiga Gosar Stran 14/49 Seminar: Numerične metode wi 1  wi  25 1408 2197 1 k1  k3  k4  k5 216 2565 4104 5 Uporabljena lokalno napako 5. reda: wi 1  wi  16 6656 28561 9 2 k1  k3  k4  k5  k6 135 12825 56430 50 55 Koeficienti k1,…,k6 so enaki: k1  hf  ti , wi  h 1   k2  hf  ti  , wi  k1  4 4   3h 3 9   k3  hf  ti  , wi  k1  k2  8 32 32   12h 1932 7200 7296   k4  hf  ti  , wi  k1  k2  k3  13 2197 2197 2197   439 3680 845   k5  hf  ti  h, wi  k1  8k2  k3  k4  216 513 4104   h 8 3544 1859 11   k6  hf  ti  , wi  k1  2k2  k3  k4  k4  2 27 2565 4104 40   Vrednost q je potrebna za določevanje primernosti začetne izbire koraka h v i-ti iteraciji in s tem pogojuje ponovni izračun z uporabo koraka qh. Potrebna je tudi za napoved začetnega koraka h pri (i+1)-ni iteraciji. 3.2 Definicija naloge (str. 288, naloga 3d) Z uporabo Runge-Kutta-Fehlberg-ove metode s toleranco TOL= 10-4, hmax=0.5 in hmin=0.05 izpiši tabelo funkcije z začetnimi pogoji problema. Dobljene rezultate metode primerjaj z eksaktnimi vrednostmi. d) y   t  2t 3.3 3 y 3  ty , 0t  2, 1 y  0  ; 3  eksaktn rešitev y  t   3  2t  6e 2 t2   1 2 Programska koda – Mathematica (*RUNGE-KUTTA-FEHLBERGOV ALGORITEM*) F[t_,y_]:=(t+2*t^3)*y^3-t*y; (*enačba funkcije*) a=0 ; (*levi del neenačbe*) b=2; (*desni del neenačbe*) =1/3; (*začetna vrednost*) tol=10^(-6); (*napaka*) hmin=0.05; (*najmanjša vrednost koraka*) hmax=0.5; (*največja vrednost koraka*) Žiga Gosar Stran 15/49 Seminar: Numerične metode (*pogoj neenačbe*) If[ab,Input["Leva vrednost mora biti nizja od desne vrednosti neenacbe\n Pritisni 1 [enter] za nadaljevanje\n"],OK=1;]; (*pogoj podajanja napake*) If[tol0,Input["Tolerancno obmocje mora biti pozitivno\n Pritisni 1 [enter] za nadaljevanje\n"],OK=1;]; (*pogoj podajanja koraka*) If[hmin<hmax&&hmin>0,OK=1,Input["Minimalna razdalja mora biti pozitivno, realno, stevilo in manjsa kot maksimalna razdalja, hmax\n Pritisni 1 [enter] za nadaljevanje\n"];]; (*izpis rezultata*) If[OK1,FLAG=Input["Izberi izpis rezultatov\n 1. Ekran\n 2. Datoteka\n Enter 1 or 2\n"]; If[FLAG2,NAME=InputString["Podaj ime datoteke\n Primer: output.dta\n"]; OUP=OpenWrite[NAME,FormatTypeOutputForm],OUP="stdout";]; Write[OUP,"Runge-Kutta-Fehlbergova metoda"]; Write[OUP,"t(i) w(i) h R"]; (*Korak 1*) H=hmax; t=a; W=; Write[OUP,t," OK=1; ",N[W,9]," 0 0 "]; (*Korak 2 in 3*) While[t<b&&OK1,K1=H*F[t,W]; K2=H*F[t+H/4,W+K1/4]; K3=H*F[t+3*H/8,W+(3*K1+9*K2)/32]; K4=H*F[t+12*H/13,W+(1932*K1-7200*K2+7296*K3)/2197]; K5=H*F[t+H,W+439*K1/216-8*K2+3680*K3/513-845*K4/4104]; K6=H*F[t+H/2,W-8*K1/27+2*K2-3544*K3/2565+1859*K4/4104-11*K5/40]; (*Korak 4*) R=Abs[K1/360-128*K3/4275-2197*K4/75240.0+K5/50+2*K6/55]/H; (*Korak 5*) If[Rtol, (*Korak 6-Aproksimacija*) t=t+H; W=W+25*K1/216+1408*K3/2565+2197*K4/4104-K5/5; (*Korak 7*) Write[OUP,N[t,9]," ",N[W,9]," ",N[H,9]," ",N[R,9]];]; (*Korak 8*) If[R>1.0 10^-20,=.84*Exp[0.25*Log[tol/R]],=10.0;]; (*Korak 9-Racunaj novi H*) If[0.1,H=.1*H,If[4,H=4.0*H,H=*H;];]; (*Korak 10*) If[H>hmax,H=hmax;]; (*Korak 11*) If[H<hmin,OK=0,If[t+H>b,If[Abs[b-t]<tol,t=b,H=b-t;];];];]; Žiga Gosar Stran 16/49 Seminar: Numerične metode If[OK0,Write[OUP,"Minimalni korak (h) prevelik\n"];]; (*Korak 12-Proces je končan*) If[OUP"OutputStream[",NAME," 3]",Print["Datoteka: ",NAME," narejena uspesno\n"]; Close[OUP]];]; 3.3.1.1 Rešitev Runge-Kutta-Fehlbergova metoda "t(i) w(i) h R" 0 0.398605 0.333333333 0.31082 0 0.398605 0 4.24957×10-7 0.683726 0.272057 0.285121 4.87189×10-7 0.970397 0.222119 0.286671 4.58196×10-7 1.26308 0.167198 0.292685 4.26611×10-7 1.56729 0.113308 0.304209 5.27336×10-7 1.73331 0.0876899 0.166017 3.90042×10-7 1.90977 0.0644823 0.176463 9.89429×10-7 2. 0.0543454 0.0902302 9.43252×10-8 Žiga Gosar Stran 17/49 Seminar: Numerične metode 4 Naloga 4 4.1 Numerično reševanje diferencialnih enačb  , a t b Diferencialne enačbe višjega reda y m t   f t , y, y... y m1 z začetnimi pogoji y  a   1 , y  a    2 ,..., y enačb prvega reda. m1 Naj bo u1  t   y  t  , u2  t   y  t  ,..., um  t   y  a   m m1 (4.1) lahko prevedemo na sistem diferencialnih t  , kar privede do sistema diferencialnih enačb prvega reda, in sicer: du1 dy   u2 dt dt du2 dy   u3 dt dt dum1 dy m2   um dt dt in   dum dy m1 m m1   y   f t , y, y... y   f  t , u1 , u2 , dt dt , um  z začetnimi pogoji u1  a   y  a   1 , u2  a   y  a    2 ,..., um  a   y m1  a    m Sedaj lahko uporabimo metodo Runge-Kutta za reševanje m-tih diferencialnih enačb z obliko du1  f1  t , u1 , u2 , dt du2  f2  t , u1 , u2 , dt dum  fm  t , u1 , u2 , dt , um  , um  , um  za a  t  b z začetnimi pogoji definiranih v izrazih (4.1). Cilj te metode je poiskati m funkcij u1,u2,…,um, katere zadovoljijo vsako diferencialno enačbo skupaj z vsemi začetnimi pogoji. Razdelimo interval [a,b] na N podintervalov z mrežo točk tj= a+ jh za vsak j= 0,1…N in h= (b-a)/N. Naj bo wi,j aproksimacija funkcije ui(tj) za vsak j=0,1,…,N in i=1,2,…,m. Definirajmo začetne pogoje kot Žiga Gosar Stran 18/49 Seminar: Numerične metode w1,0  1 , w2,0   2 ,..., wm,0   m k1,i  h fi  t j , w1, j , w2, j ,..., wm, j  za vsak i  1, 2,..., m Aproksimacijo funkcij ui preko vseh mrežnih točk rešujemo po naslednjem postopku: 1 1 1 1   k2,i  h fi  t j  , w1, j  k1,1 , w2, j  k1,2 ,..., wm, j  k1,m  za vsak i  1, 2,..., m 2 2 2 2   1 1 1 1   k3,i  h fi  t j  , w1, j  k2,1 , w2, j  k2,2 ,..., wm, j  k2,m  za vsak i  1, 2,..., m 2 2 2 2   k4,i  h fi  t j  h, w1, j  k3,1 , w2, j  k3,2 ,..., wm, j  k1,m  za vsak i  1, 2,..., m in še wi , j 1  wi , j  1  k1,i  2k2,i  2k3,i  k4,i  za vsak i  1, 2,..., m 6 Pri tem poteka j=0,1,…,N ter pazimo, da za računanje wi,j+ 1 že poznamo v prejšnji točki aproksimirano vrednost wi,j. 4.2 Definicija naloge (str 323, naloga 6) Nihalo na sliki 4-1 je dolžine 2 čevlja. Pri tem je težnostni pospešek g = 9.81 m/s2. S korakom h = 0,1 sekunde, primerjajte kot θ za naslednja začetna problema: a) b) d 2 g  sin   0 , dt 2 L d 2 g    0, dt 2 L   0    0   6  6 , ,    0  0 ,    0  0 , pri času t = 0, 1 in 2 s. 4.3 Programska koda – Mathematica a) del naloge (* Sistem dveh diferencialnih enacb *) f1[w2_]:=w2; (*w2 - kotna hitrost*) f2[t_,w1_,w2_]:=-g/L*Sin[w1]; α1=N[Pi/6]; α2=0; g=9.81; L=2; a=0;b=2;h=0.1; n=(b-a)/h; Žiga Gosar Stran 19/49 Seminar: Numerične metode t=a; (* zacetni pogoji *) w1=α1; w2=α2; Print[Row[{t,w1,w2},","]] Do[ { w1=w1, w2=w2, K11=h f1[w2], K12=h f2[t,w1,w2], K21=h f1[w2+ K12/2], K22=h f2[t+h/2,w1+K11/2,w2+K12/2], K31=h f1[w2+ K22/2], K32=h f2[t+h/2,w1+K21/2,w2+K22/2], K41=h f1[w2+ K32], K42=h f2[t+h,w1+K31,w2+K32], w1=w1+(K11+2 K21+2 K31+K41)/6, w2=w2+(K12+2 K22+2 K32+K42)/6, t=a+i h; Print[Row[{t,w1,w2},","]]},{i,1,n}] b) del naloge Ponovimo postopek iz a) dela naloge, le da upoštevamo: f2[t_,w1_,w2_]:=-g/L*w1 4.3.1.1 Rešitev Tabela 6: Rešitev z metodo Runge-Kutta Čas 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 Žiga Gosar Primer a Kot (sin()) [°] kotna hitrost 0.5235990 0.0000000 0.5113800 -0.2435090 0.4752440 -0.4766050 0.4167520 -0.6889350 0.3384870 -0.8704270 0.2440090 -1.0117700 0.1377380 -1.1051600 0.0247662 -1.1451100 -0.0894133 -1.1292000 -0.1992350 -1.0584000 -0.2993980 -0.9369440 -0.3851600 -0.7717320 -0.4525730 -0.5715260 Primer b Kot () [°] kotna hitrost napaka upoštevanja 0.5235990 0.0000000 0.00000 0.5108100 -0.2547260 0.00057 0.4730690 -0.4970080 0.00217 0.4122200 -0.7150120 0.00453 0.3312350 -0.8980890 0.00725 0.2340700 -1.0373000 0.00994 0.1254710 -1.1258300 0.01227 0.0107437 -1.1593800 0.01402 -0.1045090 -1.1362800 0.01510 -0.2146560 -1.0576900 0.01542 -0.9274280 0.01492 -0.3143170 -0.3986240 -0.7518630 0.01346 -0.4634600 -0.5395720 0.01089 Stran 20/49 Seminar: Numerične metode 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 -0.4986230 -0.5213040 -0.5196450 -0.4937160 -0.4446300 -0.3745310 -0.2865660 -0.1848020 -0.3461450 -0.1058610 0.1389240 0.3777950 0.6003590 0.7963790 0.9561400 1.0710800 -0.5056560 -0.5231520 -0.5150930 -0.4818730 -0.4251140 -0.3475900 -0.2530870 -0.1462210 -0.3009250 -0.0475782 0.2080920 0.4535970 0.6769450 0.8672240 1.0151400 1.1134700 0.00703 0.00185 -0.00455 -0.01184 -0.01952 -0.02694 -0.03348 -0.03858 Napaka 0.02 0.01 napaka 0.00 -0.01 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 -0.02 -0.03 -0.04 -0.05 čas [s] Diagram 2: Razlika med lineariziranim in nelineariziranim kotom Žiga Gosar Stran 21/49 Seminar: Numerične metode Nihalo - upoštevanje kota 0.60 0.40 kot [°] 0.20 0.00 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 -0.20 -0.40 -0.60 čas [s] RK sin(theta) RK (theta) Diagram 3: Razlika med linearizirano in nelinearizirano diferencialno enačbo Nihalo - kotna hitrost 1.5 1.0 kotna hitrost 0.5 0.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 -0.5 -1.0 -1.5 čas [s] RK sin(theta) RK (theta) Diagram 4: Razlika kotnih hitrosti med linearizirano in nelinearizirano diferencialno enačbo Žiga Gosar Stran 22/49 Seminar: Numerične metode 5 Naloga 5 Teorija aproksimacij in Čebišev polinom 5.1 Množica ortogonalnih polinomov na [-1,1] z utežno funkcijo   x  polinomi1. Lahko jo definiramo s pomočjo eksplicitne formule Tn  x  cos  n arccos x , 1 1  x2 se imenuje Čebiševi n  0,1,..., N Upoštevamo enačbo: cos  n  1   cos  n  1   2cos  cos n , n Sklepamo, da so Tn polinomi stopnje n. Z upoštevanjem superpozicije   arccos x , dobimo  T  x T  x 1 1 i j dx 1  x2   cos i cos j d  0  1   cos i cos j d   ij .  2  2 Od tu nepostredno sledi ortogonalnost Čebiševih polinomov. Slika 5-1: Čebiševi polinomi T5 in T12 Čebišev polinom ima sledeče lastnosti:  n  0,1,...  Čebišev polinom Tn na intervalu [-1,1] n različnih ničtih točk  Čebišev polinom Tn na intervalu [-1,1] ima n+1 različnih točk  1 Tn  x  1 , x  1,1 ,  2k  1  xk  cos    , k  1,..., n  2n   k  cos k , n k  0,1,..., n T0  x  1 , T1  x  x za n  1, 2,... velja rekurzivna funkcija v katerih ima globalne maksimume in minimume Oznaka polinoma »T« pride iz angleškega prevoda priimka Čebišev: Tchebycheff Žiga Gosar Stran 23/49 Tn1  x  2 xTn  x  Tn1  x , 5.2 Seminar: Numerične metode n  1, 2,... Definicija naloge (str. 516, naloga 1c) Uporabi ničle T3 za izdelavo interpolacijskega polinoma druge stopnje za dano funkcijo na intervalu [-1,1]. f  x  ln  x  2  5.3 Programska koda – Mathematica (*ČBEIŠEV POLINOM*) f[x_]=Log[x+2]; n=2; (*stopnja polinoma*) tk_=Cos[((2 n +1-2 k))/(2 n +2)]; (*ničle*) XY=Table[N[{tk,f[tk]}],{k,0,n}]; (*izračun ničel*) P2[x_]=Expand[InterpolatingPolynomial[XY,x]]; Needs["Graphics`Colors`"]; dots=ListPlot[XY,PlotStyle{Red,PointSize[0.02]},DisplayFunctionIdentity]; gr1=Plot[{f[x],P2[x]},{x,- 1.5,1.5}, PlotStyle{Red,Blue}, DisplayFunctionIdentity]; Show[gr1,dots,PlotRange{{-1.0,1.0},{0.0,1.2}}, AxesLabel{x,y},DisplayFunction$DisplayFunction] Print["f[x] = ",f[x]]; Print["Čebiševe ničle = ", XY]; Print["Interval aproksimacije je v mejah [-1.0,1.0]."]; Print["P2[x] = ", P2[x]]; (*Napaka Čebiševeka interpolacijskega polinoma P2*) e2[x_]=f[x]-P2[x]; Plot[e2[x],{x,-1.0,1.0},PlotStyle{Blue},AxesLabel{x,y}] (*Print["f[x] = ",f[x]];*) (*Print["P2[x] = ",P2[x]];*) Print["Interval aproksimacije je v mejah [-1.0,1.0]."]; Print["Graf napake e2[x]=f[x]-P2[x]"]; extrema1=FindMinimum[e2[x],{x,0.5}]; min1=extrema1[[1]]; extrema2=-FindMinimum[-e2[x],{x,-0.5}]; max1=extrema2[[1]]; extrema={max1,min1}; maxerr2=Max[Abs[{max1,min1}]]; Print["Ekstrema za e2[x]; ", extrema]; Print["|e2[x]|  ", maxerr2]; Žiga Gosar Stran 24/49 Seminar: Numerične metode 5.3.1.1 Rešitev y 1.2 1.0 0.8 0.6 0.4 0.2 x 1.0 f[x] = 0.5 0.5 1.0 Log[2+x] Čebiševe točke = {{-0.866025,0.125729},{0.,0.693147},{0.866025,1.05293}} Interval aproksimacije je v mejah [-1.0,1.0]. P2[x] = 0.693147 + 0.535318 x - 0.138426 x2 y 0.015 0.010 0.005 x 1.0 0.5 0.5 1.0 0.005 0.010 0.015 0.020 Interval aproksimacije je v mejah [-1.0,1.0]. Graf napake e2[x]=f[x]-P2[x] Ekstrema za e2[x]; |e2[x]|  Žiga Gosar {0.0147259,-0.00995048} 0.0147259 Stran 25/49 Seminar: Numerične metode Tabela 7: Aproksimativne vrednosti Chebishevega polinoma x -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 P2[x] = f(x) 0.00000000 0.09531018 0.18232156 0.26236426 0.33647224 0.40546511 0.47000363 0.53062825 0.58778666 0.64185389 0.69314718 0.74193734 0.78845736 0.83290912 0.87546874 0.91629073 0.95551145 0.99325177 1.02961942 1.06471074 1.09861229 f(x)aproks 0.01940327 0.09923603 0.17630025 0.25059596 0.32212314 0.39088179 0.45687192 0.52009352 0.58054660 0.63823115 0.69314718 0.74529468 0.79467366 0.84128412 0.88512605 0.92619945 0.96450433 1.00004069 1.03280852 1.06280782 1.09003860 f(x) - f(x)aproks 0.01940327 0.00392585 0.00602130 0.01176831 0.01434910 0.01472590 0.01313171 0.01053473 0.00724007 0.00362273 0.00000000 0.00335734 0.00621630 0.00837499 0.00965731 0.00990872 0.00899289 0.00678891 0.00318910 0.00190292 0.00857369 0.693147 + 0.535318 x - 0.138426 x2 Interpolacijski polinom 2. stopnje (rešitev v literaturi [1]): P2  x  1.052926  0.4154370  x  0.8660254  0.1384262x  x  0.8660254 P2  x  0.693147  0.535318x  0.138426 x2 Žiga Gosar Stran 26/49 Seminar: Numerične metode 6 Naloga 6 6.1 Teorija paproksimacij in Fast Foutier Transformacija Vrednosti obravnavanih funkcij lahko aproksimiramo tudi z uporabo sinusnih in kosinusnih funkcij z uporabo takoimenovane Fourier-ove vrste. Aproksimacijo k diskretno podanim točkam  x , y  2 m1 j j j 0 Sn ( x)  funkcij f  xj   y j , kjer so xj     j m  in je f  C   ,   , zapišemo v obliki: n 1 a0  a n cos  nx    a k cos  kx  bk sin  kx  2 k 1 y  Sn  x j  Koeficiente ak in bk, ki minimizirajo napako: E (a o ,..., a n , b1 ,..., bn1 )  2 m1 j 0 j  2 Zapišemo kot: 1 2 m1  y j cos  kxj  , m j 0 ak  1 2 m1  y j sin  kxj  , m j 0 bk  k  0,1, 2,..., n (6.1) k  1, 2,..., n  1 (6.2) Limito Sn(x), ko gre n   imenujemo Fourier-jeva vrsta funkcije f. V primeru, ko izvajamo interpolacijo skozi diskretno podane točke zapišemo interpolacijski polinom kot: Sm  x   x , y  a 0  a m cos  mx m1    a k cos  kx  bk sin  kx  2 k 1 2 m1 j j j 0 funkcije f, (6.3) Koeficienti a k in bk pa so enaki kot (6.1) in (6.2). Računanje trigonometrijske interpolacije po enačbah (6.2) in (6.3) imenujemo direktno reševanje, pri čemer je potrebno izvesti približno (2m)2 množenj in (2m)2 seštevanj. Druga metoda izračuna trigonometrijskega interpolacijskega polinoma potrebuje ob ustrezni izbiri m le O(mlog2m) množenj in O(mlog2m) seštevanj. To metodo imenujemo Cooley-Tukey metoda ali hitra Fourier-jeva transformacija. Pri tej metodi namesto računanja konstant a k in bk računamo kompleksne koeficiente v: 1 2 m1 ikx  ck e , m j 0 Kjer so: ck  ye 2 m1 j 0 j ik j m k  1, 2,..., n  1 , k  0,1,..., 2m 1 Iz konstant ck izračunamo koeficiente a k in bk s pomočjo Euler-jeve relacije: eiz  cos  z  i sin  z  , Žiga Gosar k  0,1,..., m Stran 27/49 Seminar: Numerične metode Za vsak k = 0,1,…,m tako izpeljemo: a k  ibk  1  k m ck . Če upoštevamo relacijo: eni   1 , n lahko za vsak n vpeljemo redukcijo potrebnih izračunov. Ker je: ck  ck  M  in je:  y e  1  e  2 m1 j 0 ik j m ij j 2, če je j sodo 1  e ij   0, če je j liho je zato: ck  ck  M  2 y2e m1 j 0 ik j  m 2  ck  ck  M  2eik m  y2 j 1e m1 j 0 (6.4) ik j  m 2  (6.5) Iz enačb (6.4) in (6.5) lahko pridobimo vse koeficiente ck. Namesto(2m)2 kompleksnih množenj je potrevno po enačbah (6.4) in (6.5) izvesti le 2m2+m množenj. Tehniko reduksije lahko vpeljemo rkrat, kjer je r = p+1 in je m = 2p. Zreducirano število kompleksnih množenj je torej enako: m2  mr 2r  2 6.2 Definicija naloge (str. 547, naloga 5c) 3c) Uporabi FFT algoritem za izračun trigonometričnega interpolacijskega polinoma 4. stopnje na intervalu [-π, π] za sledečo funkcijo: f  x  cos  x  2sin  x 5c) Uporabi aproksimacijo dobljeno v nalogi 3 za aproksimiranje sledečeka integrala in primerjaj rezultata z dejansko vrednostjo.    cos  x  2sin  xdx   6.3 Programska koda – Mathematica (*FAST FOURIER ALGORITEM*) (*računaš v intervalu [-Pi,Pi]*) Žiga Gosar Stran 28/49 Seminar: Numerične metode F[x_]:=Cos[Pi*x]-2*Sin[Pi*x]; m=4; (*stopnja trigonometričnega interpolacijskega polinoma*) n=2*m; If[n>0,For[JJ=1,JJn,JJ++,Z=-Pi+(JJ-1)*Pi/m; Y[JJ-1]=F[Z];]; OK=1,Input["Število mora biti pozitivno\n Pritisni 1 [enter] za nadaljevanje\n"];]; If[OK1,FLAG=Input["Izberi vrsto izpisa rezulta\n 1. Ekran\n 2. V datoteko\n Vnesi 1 ali 2\n"]; If[FLAG2,NAME=InputString["Podaj ime datoteke\n Priemer: output.dta\n"]; OUP=OpenWrite[NAME,FormatTypeOutputForm],OUP="stdout"]; Write[OUP,"FAST FOURIER TRANSFORMACIJA\n"]; TW=Log[2.0]; (*korak 1*)(*zamenjam: N2 za m*) N2=m; (*korak 2*) For[JJ=1,JJn,JJ++,c[JJ-1]=Y[JJ-1];]; Z=n; p=N[Round[Log[n]/TW]]; q=p-1; YY=2*Pi*I/n; =Exp[YY]; (*korak 3*) For[JJ=1,JJN2,JJ++,X=1; YY=1; For[J=1,JJJ,J++,YY=X*; X=YY;]; [JJ-1]=X; [N2+JJ-1]=-X;]; (*korak 4*) K=0; [-1]=1; (*korak 5*) For[L=1,Lp,L++, (*korak 6*) While[K<n-1, (*korak 7*) For[JJ=1,JJN2,JJ++, (*korak 8*) Z=Exp[q*TW]; M1=Round[Z]; M1=Floor[K/M1]; (*IBR - bit reversal*) R1=M1; S=0; For[RR=1,RRp,RR++,R2=Floor[R1/2]; S=2*S+R1-2*R2; R1=R2;]; Žiga Gosar Stran 29/49 Seminar: Numerične metode NP=S; (*T1 je eta*) T1=c[K+N2]; (*korak 9*) If[NP0,X=T1; T1=X*[NP-1];]; c[K+N2]=c[K]-T1; c[K]=c[K]+T1; (*korak 10*) K=K+1;]; (*korak 12*) K=K+N2;]; K=0; N2=Floor[N2/2]; q=q-1;]; (*korak 13*) While[K<n-1, (*korak 14*) R1=K; JJ=0; For[RR=1,RRp,RR++,R2=Floor[R1/2]; JJ=2*JJ+R1-2*R2; R1=R2;]; (*korak 15*) If[JJ>K,T3=c[K]; c[K]=c[JJ]; c[JJ]=T3;]; (*korak 16*) K=K+1;]; (*koraka 17 in 18*) Write[OUP,"Koeficienti c(0), ... , c(2m-1)"]; For[JJ=1,JJn,JJ++, YY=-(JJ-1)*Pi*I; X=Exp[YY]; YY=X*c[JJ-1]; c[JJ-1]=N[YY/(0.5*n)]; K=JJ-1; Write[OUP,K," ",c[JJ-1]];]; Write[OUP,"\n"]; Write[OUP,"Koeficienti a(0), ... , a(m)"]; For[JJ=1,JJm+1,JJ++, Write[OUP,N[Re[c[JJ-1]],9]];]; Write[OUP,"\n"]; Write[OUP,"Koeficienti b(1), ... , b(m-1)"]; For[JJ=2,JJm,JJ++, Write[OUP,N[Im[c[JJ-1]],9]];]; (*korak 19*) If[OUP"Izhod[",NAME," 3]",Print["Izhodna datoteka: ",NAME," je uspešno narejena.\n"]; Close[OUP];];]; 6.3.1.1 Rešitev: FAST FOURIER TRANSFORM Žiga Gosar Stran 30/49 Seminar: Numerične metode Koeficienti c(0), ... , c(2m-1) 0 1 2 -0.126426 0.260272 -0.102219  -0.301114+0.275406  3 1.12137 -2.05295  4 0.0458965 5 6 7  1.12137 +2.05295  -0.301114-0.275406  0.260272 +0.102219  Koeficienti a(0), ... , a(m) -0.126426 0.260272 -0.301114 1.12137 0.0458965 Koeficienti b(1), ... , b(m-1) -0.102219 0.275406 -2.05295 6.3.1.2 Rešitev 3c) [1] P  x3c  0.1264264  0.2602724  cos  x   0.3011140  cos  2  x   1.121372  cos  3  x   0.04589648  cos  4  x   0.102219  sin  x   0.2754062  sin  2  x   2.052955  sin  3  x  3c) izračunano s pomočjo spisanega programa P  x5c  0.126426  0.260272  cos  x   0.301114  cos  2  x  1.12137  cos  3  x   0.0458965  cos  4  x   0.102219  sin  x   0.275406  sin  2  x   2.05295  sin  3  x  Žiga Gosar Stran 31/49 Seminar: Numerične metode FFT 2.5 1.5 0.5 -2.2 -1.2 -0.2 -0.5 y -3.2 0.8 1.8 2.8 -1.5 -2.5 -3.5 x - interval od -PI do PI Naloga 3c Naloga 5c napaka Naloga 5c) Rešitev naloge 3c) NIntegrate[-0.1264264+0.26027239 Cos[x]-0.301113969 Cos[2 x]+1.1213715 Cos[3 x]+0.045896473 Cos[4 x]-0.10221901 Sin[x]+0.27540621 Sin[2 x]-2.052954949 Sin[3 x],{x,-Pi,Pi}] -0.79436 Eksaktna vrednost integrala   cos  x  2sin  xdx   NIntegrate[Cos[Pi x]-2 Sin[Pi x], {x,-Pi,Pi}] -0.273938 Rešitev naloge 5c se povsem ujema z rešitvijo podano v [1] na strani 811. Žiga Gosar Stran 32/49 Seminar: Numerične metode 7 Naloga 7 7.1 QR algoritem Za računanje lastnih vrednosti simetrične realne matrike A reda n x n je zelo znan QR algoritem. Matriko A je mogoče razcepiti na produkt ortogonalne matrike Q in zgornje trikotne matrike R, A=QR. Pomnožimo enačbo z leve s Q-1=Q+t in z desne s Q. Dobimo QtAQ=RQ. Torej sta matriki podobni, zato imata enake lastne vrednosti. Na tej lastnosti temelji QR algoritem za računanje lastnih vrednosti matrike [2]. Naj bo A realna simetrična matrika, označimo jo z A1. Če želimo izračunati njene lastne vrednosti, jo razcepimo v produkt Q1R1. Nato izračunamo A2=R1Q1. Zopet razcepimo matriko A2 v produkt Q2R2 in izračunamo A3=R2Q2 itd. 7.1.1 QR iteracija S QR iteracijo na ekonomičen način izračunamo vse lastne vrednosti matrike. Za redukcijo na tridiagonalno obliko zaradi simetrije porabimo približno polovico toliko operacij kot za nesimetrično matriko, a to še vedno pomeni O(n3) operacij. Med samo QR iteracijo je razlika večja. Zaradi tridiagonalne oblike lahko en korak QR iteracije izvedemo z zahtevnostjo O(n). Pri QR iteraciji najprej poiščemo ortogonalno matriko Q, da je T = QAQT tridiagonalna, potem pa delamo QR z enojnim premikom: T0 = T k = 0, 1, . . . izberi premik σk Tk − σk I = QkRk (izračunaj QR razcep) Tk+1 = RkQk + σk I Naj bo  a1 k   k  b1 Tk      b1 k  a 2 k  b2 k  bn 2 k a n 1 k bn 1 k      k bn 1   k  an   a  k  b k   • Wilkinsonov premik oz. dvojni premik: za σk vzamemo tisto lastno vrednost matrike  nk1 n k1  ,      bn 1 a n  ki je bližja a n k  . Kako izberemo premik:   Obdiagonalne elemente postavimo na 0, kadar velja bj k    a jk   a jk1 Žiga Gosar Stran 33/49 Seminar: Numerične metode 7.2 Definicija naloge (str. 595, naloga 2d) Z uporabo QR algoritma določi, v tolerančnem območju 10-5, vse lastne vrednosti sledeče matrike. 0 0  5 1  1 4.5 0.2 0   0 0.2 1 0.4  0 0.4 3 0  0 0 0 1 7.3 0 0  0  1 3 Programska koda – Mathematica (*QR ALGORITEM*) A={{5,-1,0,0,0},{-1,4.5,0.2,0,0},{0,0.2,1,-0.4,0},{0,0,-0.4,3,1},{0,0,0,1,3}}; (*Print["Podana je matrika"];*) (*Print["A = ",MatrixForm[A]];*) QRmetoda[A_,m0_]:=Module[{A0=N[A],A1,i,m=m0},Print["A"0," = ", MatrixForm[Chop[A0,1.0 10-6]]] For[i=1,i<m,i++, {Q0,R0}= QRDecomposition[A0]; A1=R0.Transpose[Q0]; Print["A"i+1," = ", MatrixForm[Chop[A1,1.0 10-6]]]; A0=A1]; Return[A0];]; QRpremik[A_,_]:=Module[{B=N[A],A1=N[A],i,m=m0}, B=A1; n=Length[A1]; =Table[A1[[i,i]],{i,n}]; sum=0; For[m=n,2m,m--, B=A1[[Range[1,m],Range[1,m]]]; i=1; While[And[i10,Length[B]m,Abs[B[[m,m-1]]]>], M=B[[{m-1,m},{m-1,m}]]; =Eigenvalues[M]; j=Position[Abs[- [B[[m,m]]],Min[Abs[-B[[m,m]]]]][[1,1]]; s=[[j]]; (*s=Wilkinsonov premik*) sum=sum+s; {Q1,R1}=QRDecomposition[B-s IdentityMatrix[m]]; Print["\n","Premik ", "s"i, " = ", s]; T ", MatrixForm Transpose Chop Q1 Print["Q"i, " Print["R"i, " = ",MatrixForm[Chop[R1]]]; B=R1.Transpose[Q1]; ]; T Print["R"i,"Q"i , " = ", MatrixForm[Chop[B]]]; i=i+1;]; A1[[Range[1,m-1],Range[1,m-1]]]=B[[Range[1,m-1],Range[1,m-1]]]; A1[[m,m-1]]=A1[[m-1,m]]=0; A1[[m,m]]=sum; Print["A1 = ",MatrixForm[Chop[A1]]];]; A1[[1,1]]=sum-s+[[3-j]]; Print["\n", "Končni odgovor je:"]; Print["d = ", MatrixForm[Chop[A1]]]; Return[A1];]; d=QRpremik[A,1.0 10-6]; Žiga Gosar Stran 34/49 Seminar: Numerične metode (*g=QRmetoda[A,203];*) Print["A = ",MatrixForm[A]]; Print["{i} = ", Eigenvalues[N[A]]]; Print[""]; Print["D = ",MatrixForm[Chop[d,1.0 10-6]]]; Print["{i} = ", Eigenvalues[N[d]]]; 7.3.1.1 Rešitev Končni odgovor je: D ={{5.784, 0, 0, 0, 0}, {0, 0.8903, 0, 0, 0}, {0, 0, 3.72756, 0, 0}, {0, 0, 0, 2.07069, 0}, {0, 0, 0, 0, 4.02742}} {i} = {5.784,4.02742,3.72756,2.07069,0.8903} A203={{5.784, 0, 0, 0, 0}, {0, 4.02743, 0, 0, 0}, {0, 0, 3.72756, 0, 0}, {0, 0, 0, 2.07071, 0}, {0, 0, 0, 0, 0.8903}} Eksaktna rešitev: G={{5, -1, 0, 0, 0},{-1, 4.5, 0.2, 0, 0},{0, 0.2, 1, -0.4, 0},{0, 0, -0.4, 3, 1},{0, 0, 0, 1, 3}}; N[Eigenvalues[G]] {5.784,4.02743,3.72756,2.07071,0.8903} Žiga Gosar Stran 35/49 Seminar: Numerične metode 8 Naloga 8 8.1 Kontinuitetna metoda S kontinuitetno metodo za nelinearni sistem enačb je možno rešiti zbirko težav. Posebej za reševanje problema v obliki F(x)=0 ki ima neznano rešitev x*, predpostavimo, da je družina opisane težave z uporabo parametra λ, ki prevzamejo vrednost v intervalu [0,1]. Problem z znano rešitvijo x(0) ustreza λ=0 in je problem z neznano rešitvijo x(1)=x*, ki ustreza λ=1. Predpostavimo da je x(0) začetna aproksimacija rešitve F(x*)=0. G:0,1  Definiramo z n  n G   , x    F  x   1    F  x  -F  x  0    F  x      1 F  x  0   (8.1) Določili bomo različne vrednosti λ, kot rešitev enačbe G  , x   0 0  G   , x   F  x  -F  x  0   Ko je λ=0, enačba preide v obliko kjer je x(0) rešitev. Ko je λ=1 enačba preide v obliko 0  G 1, x   F  x  kjer je rešitev x(1)=x*. Funkcija G s parametrom λ nam zagotavlja družino funkcij, ki lahko vodijo od znane vrednosti x(0) do rešitve x(1)=x*. Funkcija G se imenuje homotopija med funkcijo G(0,x)=F(x)-F(x(0)) in funkcijo G(1,x)=F(x). Cilj kontinuitetne metode je, določiti postopek reševanja, iz znane rešitve x(0) funkcije G(0,x) do neznane rešitve x(1)=x* X1 funkcije G(1,x)=0, ki reši F(x)=0. Najprej predpostavimo, da je x(λ) edinstvene rešitve enačbe G   , x  0 , Na set   0,1 x    0   1 x    je mogoče gledati kot na krivuljo v (8.2) n od x(0) do x(1)=x*, ki je parameteriziran z λ. Kontinuitetna metoda poišče zaporedje korakov vzdolž krivulje, ki so odvisne od m k 0 kjer je λ0 = 0 < λ1 < …< λm = 1. Če sta funkciji   x    in G odvedljivi, potem odvajanje enačbe (8.2) v odvisnosti od λ k Žiga Gosar Stran 36/49 0 G   , x  λ   λ Seminar: Numerične metode  G   , x  λ   x x  λ  . Enačbo rešimo za x    , ki je sledeče oblike  G   , x  λ    G   , x  λ   . x        x λ     1 To je sistem diferencialnih enačb z začetnim pogojem x(0). G   , x  λ    F(x( ))+   1 F(x( )) Ker je lahko določlimo oboje, Jakobijevo matriko f1  f1  x  x  λ   x  x  λ   2  1 f2  f2 x  λ    x  λ  G  x2  , x  λ     x1  x    fn  x  λ   fn  x  λ    x1 x2 in G   , x  λ   λ f1  x  λ    xn  f2  x  λ    xn   J  x  λ    fn x  λ     xn  F(x(0)) . Zato sistem diferencialnih enčab lahko zapišemo kot x       J  x  λ   F  x  λ   , 1 0  λ  1, z začetnim pogojem x(0). Na splošno ima sistem diferencialnih enačb, ki jih moramo rešiti za naš kontinuitetni problem, obliko  f1  x  0    1   , x1 ,..., xn         , x ,..., x f x 0       1    2 1 2 n     J  x ,..., x  1 n  .        fn  x  0    n   , x1 ,..., xn     Da za reševanje sistema enačb uporabimo Runge-Kutta metodo 4.reda, izberemo N>0 in h=(1-0)/N. Interval [0,1] razdelimo v N podintervalov z mrežami točk λj=jh, j=0,1,…,N. Uporabimo zapis wi,j za vsak j=0,1,…,N in i=1,…,n za označevanje približka xi(λj). Začetne pogoje zapišemo w1,0= x1(0), Žiga Gosar w2,0= x2(0), … , wn,0= xn(0). Stran 37/49 Seminar: Numerične metode Predpostavimo, da smo w1,j, w2,j, … , wn,j, izračunali. w1,j+ 1, w2,j+ 1, … , wn,j+ 1, dobimo s pomočjo enačb k1,i  hi   j , w1, j , w2, j ,..., wn, j  , i  1, 2,..., n ; h 1 1 1   k2,i  hi   j  , w1, j  k1,1 , w2, j  k1,2 ,..., wn, j  k1,n  , 2 2 2 2   h 1 1 1   k3,i  hi   j  , w1, j  k2,1 , w2, j  k2,2 ,..., wn , j  k2,n  , 2 2 2 2   k4,i  hi   j  h, w1, j  k3,1 , w2, j  k3,2 ,..., wn, j  k3,n  , in končno wi , j 1  wi , j  8.2 1  k1,i  2k2,i  2k3,i  k4,i  , 6 i  1, 2,..., n ; i  1, 2,..., n ; i  1, 2,..., n ; i  1, 2,..., n . Definicija naloge (str. 642, naloga 2) Nelinearni sistem enačb f1  x1 , x2   x12  x22  2 x2  0 f2  x1 , x2   2 x1  x22  6  0 Ima dve rešitvi, (0.625204094, 2.179355825)t in (2.10951192, -1.334532188)t. Uporabi metodo Runge-Kutta četrtega reda z N=1 in aproksimiraj rezultat z: i. x  0    0, 0  t 8.3 ii. x  0   1,1 t iii. x  0    3, 2  t Programska koda – Mathematica (*KONTINUITETNA METODA ZA SISTEM NELINEARNIH ENAČB*) (*enačba 1 =x1^2-x2^2+2 x2; enačba 2 =2 x1+x2^2-6*) n=2; (*število enačb*) N1=4;(*število korakov*) For[i=1,in,i++, TEMP[i]=Input["To je Kontinuitetna metoda za sistem nelinearnih enačb.\n Če imaš več kot 2 enačbi, program dopolni.\n Podaj funkcijo F["<>ToString[i]<>"] v odvisnosti od x1,x2\n"]; F[i][x1_,x2_(*,x3_*)]:=Evaluate[TEMP[i]];]; (*dodajaj, če imaš več spremenljivk oz. enačb*) For[i=1,in,i++, DER=D[TEMP[i],x1];(*odvod enačbe po x1*) P[i,1][x1_,x2_(*,x3_*)]:=Evaluate[DER];(*dodajaj, če imaš več spremenljivk oz. enačb*) DER=D[TEMP[i],x2]; P[i,2][x1_,x2_]:=Evaluate[DER];]; If[N10,Input["Število korakov mora biti pozitivno in celo število.\n Pritisni Enter 1 [enter] za nadaljevanje."],OK=1;]; Žiga Gosar Stran 38/49 Seminar: Numerične metode For[i=1,in,i++,X[i-1]=Input["Podaj začetno vrednost X("<>ToString[i]<>")"];]; If[OK1,FLAG=Input["Izberi izpis podatkov:\n 1. Ekran\n 2. Datoteka\n Vnesi 1 ali 2"]; If[FLAG2,NAME=InputString["Input the file name\n For example: output.dta\n"]; OUP=OpenWrite[NAME,FormatTypeOutputForm],OUP="stdout";]; Write[OUP,"KONTINUITETNA MATODA ZA SISTEM NELINEARNIH ENAČB"]; (*korak 1*) H=N[1.0/N1]; RMM[1]=0.5; RMM[2]=0.5; RMM[3]=1.0; RMM[4]=0.0; For[i=1,in,i++,X1[i-1]=N[X[i-1]];]; For[i=1,in,i++,B[i-1]=N[-F[i][X[0],X[1](*,X[2]*)]];(*dodajaj, če imaš več spremenljivk oz. enačb*) B[i-1]=N[H*B[i-1]];]; K=1; (*korak 2*) While[OK1&&KN1, (*koraki 3-6*) KK=1; For[i=1,in,i++, X[i-1]=N[X1[i-1]];]; While[OK1&&KK4, For[i=1,in,i++, For[J=1,Jn,J++, A[i-1,J-1]=N[P[i,J][X[0],X[1](*,X[2]*)]];]; (*Jakobijeva matrika*) A[i-1,n]=N[B[i-1]];]; (*Rešuješ linearni sistem A K (KK)=b*) q=n-1; OK=1; i=1; While[OK1&&iq, Z=Abs[A[i-1,i-1]]; IR=i; IA=i+1; For[J=IA,Jn,J++, If[Abs[A[J-1,i-1]]>Z, IR=J; Z=Abs[A[J-1,i-1]];];]; If[Z10^-20, OK=0, If[IRi, For[J=i,Jn+1,J++, c=A[i-1,J-1]; A[i-1,J-1]=A[IR-1,J-1]; A[IR-1,J-1]=c;];]; For[J=IA,Jn,J++, c=N[A[J-1,i-1]/A[i-1,i-1]]; If[Abs[c]10^-20, c=0;]; For[L=i,Ln+1,L++, A[J-1,L-1]=N[A[J-1,L-1]-c*A[i-1,L-1]];];];]; i=i+1;]; Žiga Gosar Stran 39/49 Seminar: Numerične metode If[OK1, If[Abs[A[n-1,n-1]]10^-20, OK=0, Y[n-1]=N[A[n-1,n]/A[n-1,n-1]]; For[i=1,iq,i++, J=n-i; JA=J+1; c=A[J-1,n]; For[L=JA,Ln,L++, c=N[c-A[J-1,L-1]*Y[L-1]];]; Y[J-1]=N[c/A[J-1,J-1]];];];]; If[OK0,Print["Linearni sistem je singularen\n"];]; If[OK1, For[i=1,in,i++, RK1[KK,i-1]=N[Y[i-1]]; X[i-1]=N[X1[i-1]+RMM[KK]*RK1[KK,i-1]];];]; KK=KK+1;]; (*korak 7*) If[OK1, For[i=1,in,i++, X1[i-1]=N[X1[i-1]+(RK1[1,i-1]+2*RK1[2,i-1]+2*RK1[3,i-1]+RK1[4,i-1])/6];]; Write[OUP,K]; For[i=1,in,i++, Write[OUP,N[X1[i-1],10]];]; K=K+1;];]; (*korak 8*) If[OUP"Izhod[",NAME," 3]", Print["Datoteka: ",NAME," je bila uspešna narejena.\n"]; Close[OUP]];]; 8.3.1.1 Rešitev začetna aproksimacija N=1 1 1 2 3 korak N=4 4 rešitev v literaturi [1] N=1 napaka N=4 Žiga Gosar x(0,0) 2.303988 -2.001099 0.723918 -0.236104 1.297477 -0.639325 1.742115 -1.009680 2.109852 -1.335660 3 -2.25 0.696012 0.248901 0.890148 0.914340 x(1,1) 0.597097 2.257968 0.833515 1.444486 0.738038 1.739998 0.672260 1.977161 0.624812 2.180382 0.421052 2.618421 0.176571 0.360453 0.204285 0.438039 x(3,-2) 2.109446 -1.334563 2.799749 -1.844045 2.586471 -1.681386 2.357554 -1.511586 2.109512 -1.334532 2.713110 -1.362773 0.603664 0.028210 0.603598 0.028241 Stran 40/49 Seminar: Numerične metode 9 Naloga 9 9.1 Hiperbolične parcialne diferencialne enačbe Metoda končnih razlik je aproksimativna metoda, ki temelji na interpolacijskem pristopu. Aproksimativna rešitev je tako zasnovana na končni množici diskretnih parametrov, ki aproksimirajo vrednosti primarnih spremenljivk v končno mnogih točkah izbranega območja. Diferencialne enačbe se rešuje s pomočjo pretvorbe diferencialnih operatorjev v diferenčne, torej s pomočjo enačb za aproksimacijo odvodov. Pri tem morajo biti izpolnjeni vsi robni pogoji obravnavanega problema. Za transformacijo diferencialnih zvez v diferenčne funkcijo p(x) razvrstimo v Taylorjevo vrsto okoli točke x0 , slika 9.1. Slika 9-1: Točke območja p  x0  h   p  x0   p  x0  h h2 h3  p  x0   p  x0   1! 2! 3! Naj bo h<1. Ker je tako h3<<h, lahko člene od h2 naprej zanemarimo. Dobimo h h2 p  x0  h   p  x0   p  x0   p  x0  1! 2! in p  x0  h   p  x0   p  x0  h h2   p  x0  1! 2! Zgornja izraza zapišemo kot w1  w0  hw0  h2  w0 2 (9.1) w1  w0  hw0  h2  w0 2 (9.2) in Če sedaj oba izraza seštejemo, dobimo izraz za centralne razlike za drugi odvor, in sicer w0  w1  2w0  w1 h2 (9.3) Centralne razlike za prvi odvod dobimo če zgornja izraza (9.1) in (9.2) odštejemo, Žiga Gosar Stran 41/49 Seminar: Numerične metode w0  w1  w1 2h (9.4) Desne razlike za prvi odvod enostavno izrazimo iz (9.1), kjer zanemarimo člene od h naprej in dobimo w0  w0  w1 h (9.5) Na ta način določene desne razlike lahko privedejo do velikih napak, zato jih raje določimo na 4h2  podlagi seštevanja izraza (9.1) ter w2  w0  2hw0  w0 . Sledi 2 w0  3w0  4w1  w2 2h (9.6) Za primer valovne enačbe, ki je splošno podana z 2 2 p 2  p x ,  t    x, t   0 ,  t 2 x2 0 xl , in robnimi ter začetnimi pogoji p  0, t   0 in p  x, t   f  x velja p  l , t   0 , za in k h (9.7) 0t p  x, t   g  x , za t wi , j 1  2 1   2  wi , j   2  wi 1, j  wi 1, j   wi , j 1 Pri tem je   0t, 0  x 1 i  1, 2,..., m  1 j  1, 2,... (9.8) , kjer k predstavlja časovni korak, h pa korak v smeri x, h= l/m. Najprej moramo določiti P i,1. To storimo z desnimi razlikami in uporabo začetnih pogojev, in sicer z uporabo izraza (9.5) wi ,1  wi ,0  k  g  xi  , i  1, 2,..., m  1 (9.9) ki je enostavnejši in lahko prinese končnim rezultatom večje napake, oziroma wi ,1  2 1   2  f  xi   2 2  f  x   f  x    k  g  x  , i  1, 2,..., m 1 i 1 i 1 i (9.10) ki je natančnejši, saj upošteva še člen s k2. 9.2 Definicija naloge (str 725, naloga 7) Zračni tlak p(x,t) v cevi je opisan z valovno enačbo 2 p 1 2 p  , x2 c 2 t 2 0  x  1, 0t, kjer je l dolžina cevi ter c fizikalna konstanta. Če je cev odprta, potem sta robna pogoja podana z Žiga Gosar Stran 42/49 p  0, t   p0 in p  l , t   p0 . Seminar: Numerične metode Če pa je cev zaprta na koncu, kjer je x= l, potem sta rovna pogoja naslednja p  0, t   p0 in p l, t   0 . x Naj bo c = 1, l =1 in začetna pogoja p  x,0   p0 cos 2 x in p  x, 0   0 , x 0  x  1. a) Aproksimiranje tlaka za primer odprte cevi z p0 = 0.9 pri x = 0 za t =0.5 in t = 1, z uporabo algoritma 12.4 z h= k= 0.1. b) Modificirajte Algoritem 12.4 za primer zaprte cevi z p0 = 0.9 ter aproksimirajte p(0.5,0.5) in p(0.5,1) za h= k= 0.1. 9.2.1 Programska koda – Mathematica a) del naloge -----------------------------------------------------------------(* Podatki *) c=1; l=1; p0=0.9; f[x_]:=p0*Cos[2*Pi*x]; g[x_]:=0; k=h=0.1; Tmax=2; (* Maksimalni cas *) λ=c*k/h; m=Round[l/h]; n=Round[Tmax/k]; w=Array[W,{m+1,n+1}]; (* Tocke mreze *) (*korak 3*) Do[{ w[[1,j]]=p0, w[[m+1,j]]=p0},{j,2,n+1}] w[[1,1]]=f[0]; w[[m+1,1]]=f[l]; (*korak 4*) Do[{ w[[i,1]]=f[(i-1)*h], w[[i,2]]=(1-λ^2)*f[(i-1)*h]+λ^2/2*(f[i*h]+f[(i-2)*h])+k*g[(i-1)*h]},{i,2,m}] Žiga Gosar Stran 43/49 Seminar: Numerične metode (*korak 5*) Do[{ Do[{ w[[i,j+1]]=2*(1-λ^2)*w[[i,j]]+λ^2*(w[[i+1,j]]+w[[i-1,j]])-w[[i,j-1]]}, {i,2,m}]}, {j,2,n}] (*korak 6*) Print["p(0.5,0.5)"] w[[0.5/h+1,0.5/k+1]] Print["p(0.5,1.0)"] w[[0.5/h+1,1.0/k+1]] Rezultati za a) del naloge p(0.5,0.5)=0.9 p(0.5,1.0)=2.7 Zrači tlak v odprti cevi 3.0 2.5 2.0 1.5 1.0 p 0.5 N11 0.0 N9 -0.5 N7 x N5 -1.0 N3 N1 1 2 3 4 5 6 7 8 21 18 19 20 15 16 17 12 13 14 11 10 9 t b) del naloge V tem primeru uporabimo še leve razlike za prvi odvod, s katerim računamo zračni tlak na koncu cevi (x= l) za vsak časovni korak, wm, j  4wm1  wm2  2  h  g  l  , 3 j  1, 2,... (9.11) -----------------------------------------------------------------(* Podatki *) Žiga Gosar Stran 44/49 Seminar: Numerične metode c=1; l=1; p0=0.9; f[x_]:=p0*Cos[2*Pi*x]; g[x_]:=0; k=h=0.1; Tmax=2; (* Maksimalni cas *) λ=c*k/h; m=Round[l/h]; n=Round[Tmax/k]; w=Array[W,{m+1,n+1}]; (* Tocke mreze *) (*korak 3*) Do[{w[[1,j]]=p0},{j,2,n+1}] w[[1,1]]=f[0]; w[[m+1,1]]=f[l]; (*korak 4*) Do[{ w[[i,1]]=f[(i-1)*h], w[[i,2]]=(1-λ^2)*f[(i-1)*h]+λ^2/2*(f[i*h]+f[(i-2)*h])+k*g[(i-1)*h] },{i,2,m}] w[[m+1,2]]=(4*w[[m,2]]-w[[m-1,2]]+2*h*g[m*h])/3; (*9.11*) Do[{ Do[{ (*korak 5*) w[[i,j+1]]=2*(1-λ^2)*w[[i,j]]+λ^2*(w[[i+1,j]]+w[[i-1,j]])-w[[i,j-1]] },{i,2,m}], w[[m+1,j+1]]=(4*w[[m,j+1]]-w[[m-1,j+1]]+2*h*g[m*h])/3; (*9.11*) },{j,2,n}] (*korak 6*) Print["p(0.5,0.5)"] w[[0.5/h+1,0.5/k+1]] Print["p(0.5,1.0)"] w[[0.5/h+1,1.0/k+1]] Rezultati za b) del naloge p(0.5,0.5)=0.9 p(0.5,1.0)=0.918793 Žiga Gosar Stran 45/49 Seminar: Numerične metode Zračni tlak na koncu cevi 3.0 2.5 2.0 1.5 1.0 p 0.5 0.0 -0.5 -1.0 N11 N9 N7 x N5 N3 N1 1 Žiga Gosar 3 5 7 9 11 13 15 17 19 21 t Stran 46/49 Seminar: Numerične metode 10 Primerjava numeričnih metod Odločil sem se, da bom med seboj primerjal metode s katerimi sem reševal probleme v moji seminarski nalogi. Med seboj bom primeral Taylorjev polinom 2. in 4. reda, metodo Runge-Kutta 4. reda in Runge-Kutta-Fehlbergovo metodo, ki ima spremenljivi korak. Z metodami bom rešil primer naloge z naslednjimi podatki: y  y  t 2  1 , y  0   0,5 0t  2, Kjer je h  0, 2 , N  10 , tol  105 , hmax  0, 25 in hmin  0,01 . t 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 y(t) 0.5000000 0.8292986 1.2140877 1.6489406 2.1272295 2.6408591 3.1799415 3.7324000 4.2834838 4.8151763 5.3054720 T(2) 0.5000000 0.8300000 1.2158000 1.6520760 2.1323327 2.6486459 3.1913480 3.7486446 4.3061464 4.8462986 5.3476843 T(4) 0.5000000 0.8293000 1.2140910 1.6489468 2.1272396 2.6408744 3.1799640 3.7324321 4.2835285 4.8152377 5.3055554 RK4 0.5000000 0.8292933 1.2140762 1.6489220 2.1272027 2.6408227 3.1798942 3.7323401 4.2834095 4.8150857 5.3053630 RKF 0.5000000 0.9204886 1.3964910 1.9537488 2.5864260 3.2604605 3.9520955 4.6308268 5.2574861 5.3054896 y(t(RKF)) 0.5000000 0.9204873 1.3964884 1.9537446 2.5864198 3.2604520 3.9520844 4.6308127 5.2574687 5.3054720 t(RKF) 0.0 0.2500000 0.4865522 0.7293332 0.9793332 1.2293332 1.4793332 1.7293332 1.9793332 2.0000000 h(RKF) 0.0 0.2500000 0.2365522 0.2427810 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.0206668 Tabela 8: Napaka metod Žiga Gosar T(2) T(4) RK4 RKF 0 0.0007014 0.0017123 0.0031354 0.0051032 0.0077868 0.0114065 0.0162446 0.0226626 0.0311223 0.0422123 0 0.0000014 0.0000033 0.0000062 0.0000101 0.0000153 0.0000225 0.0000321 0.0000447 0.0000614 0.0000834 0 0.0000053 0.0000115 0.0000186 0.0000268 0.0000364 0.0000473 0.0000599 0.0000743 0.0000906 0.0001090 0 0.0000013 0.0000026 0.0000042 0.0000062 0.0000085 0.0000111 0.0000141 0.0000174 0.0000176 Stran 47/49 Seminar: Numerične metode Primerjava metod 0.00012 0.00010 napaka 0.00008 0.00006 0.00004 0.00002 0.00000 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 t T(4) Žiga Gosar RK4 RKF Stran 48/49 Seminar: Numerične metode 11 Literatura [1] R. L. Burden and J. D. Faires, Numerical Analysis, seventh edition, Brooks/Cole, California, 2011. [2] J. Petrišič, Uvod v Matlab: za inženirje, Fakulteta za strojništvo, 2011. [3] J. Petrišič, Reševanje enačb, Fakulteta za strojništvo, 1996. [4] J. Petrišič, Interpolacija in osnove računalniške grafike, Fakulteta za strojništvo, 1999. [5] A. Kotar, Numerične metode: Zapiski z vaj, 2006 Žiga Gosar Stran 49/49