Csao!
Integral[sin(x) dx][0...Pi] -> - analitikusan Integral[-cos(x)][0...Pi] ->
-cos(Pi) - -cos(0) = 2 Ha elore tudjuk a fuggvenyt, erdemes analitikus
keplettel szamolni ;)
De neked nyilván a numerikus eljárás kell. Melegen tudom ajanlani a
Numerical recipes Internetrol is letoltheto verziojat www.nr.com. Kb 20 MB
PDF, de vannak benne igen hasznos dolgok. 4.fejezet: Integrálás. Figyelem,
sok hiba és nem korszerű eljárások rendszeresen találhatóak benne, ezért
inkább a www.NetLib.org -ot tudom javasolni, az komoly, de nehezebb rajta
elmenni (főleg FORTRAN).
Azért egyszerűbb dolgokra és/vagy bizonyos témakörökben a NumRec is jó.
Röviden amit erről ír (<> jelek között a megjegyzéseim):
- az I = Integral[f(x) dx][x1...xN] integrál legegyszerűbb megközelítésben
úgy kapható meg, hogy az abszcissza ment én több x helyen (megfelelő
gyakorisággal) összeadogatjuk az f(x) fgv értékét. (eg yéb: diffegyenlet
megoldás, függvény közelítés, FFT stb)
Jelölések, ha N darab, h szélességű részre osztjuk fel az integrálandó
tartomán yt, és a részek alatti terület összegét akarjuk megkapni:
fgv: __-----------___
-----___
| | | | | | |
-----+---+---+---+---+---+---+----
x1 x2 x3 .......... xN-1 xN
<-h->
- N+1 kiertekelesi pont lesz.
- a függvény értékek jelölése legyen f(xi)=> fi vagyis pl. f(x1) helyett f
1-et írunk.
- a feladatban az a kihívás, hogy kívánt pontossággal állítsuk elő I-t,
miközbe n a lehető legkevesebbszer értékeljük ki f(x)-et.
- vannak ún zárt és nyitott képletek az integrálszámításra. A zártak az
a,b hel yeken is felhasználják f értékét. Ez azonban néha nem lehetséges
(pl szingularitás), vagy nehézkes. Ekkor kellenek a nyitott formulák,
amelyek bentről extrapolálnak.
- a régi jó Trapéz, Simpson, Simpson 3/8, Bode stb szabályok által adott
képlet ek ma már túlhaladottnak számítanak, lehetőleg ne használjuk őket.
Egy általánosan jól használható formula a kibővített trapéz szabály, ami
gyorsa n változó függvényekre is jó: I = h * [0.5*f1 + f2 + f3 + ... +
fN-1 + 0.5*fN ] /Trap1/
A képlethez tartozó hibatag: O( f ' ' * (x1 - xN)^3 /N^2 )
vagyis a hiba a következőkkel arányos:
- f második deriváltja valahol a tartományon belül - nem tudható hol kell deriv
áltat venni.
- integrálási tartománytól köbösen - de ezen úgysem lehet változtatni.
- N, a felosztási régiók száma -> fordított négyzetesen. vagyis 2* annyi
régió 4* pontosabb eredményt ad <elvileg, kerekítési hibák nélkül>.
Az alábbiakban 100 iterációra néztem meg I-t különféle formulákkal.
I[Trap1] = 1.933
Ez túró. Lehet, hogy más formula kellene, pl ennek egy köbösen konvergáló válto
zata:
I = h * [5/12 *f1 + 13/12 * f2 + f3 + ... + 13/12* fN-1 + 5/12 *fN ] /Trap2/
Ez már O(1/N^3), és tényleg gyorsabb is.
I[Trap2] = 1.99999997
De talán van még jobb O(1/N^4): I = h * [1/3 *f1 + 4/3*f2 + 2/3*f3 +
4/3*f4 + ... + 4/3* fN-1 + 1/3*fN ] /Simp Ext1/ SimpExt1 = "Simpson
továbbfejlesztett formulája". I[SimpExt1] = 2.00000001
tehát
int iter, N = 100;
double h = (b-a)/N;
for(iter=1; iter<=N+1; iter++)
{
double mult;
if( iter==1 || iter==N+1 ) //elso, utolso
mult = 1.0/3.0;
else
if( iter%2 ) //3,5,7...
mult = 2.0/3.0;
else //2,4,6...
mult = 4.0/3.0;
I += mult*IntegFunc(a +(iter-1)*h);
}
I *= h;
és
double IntegFunc(double x)
{
return sin(x);
}
A NumRec ben még van egy kis intelligens rutin is, ami folyamatosan növeli N-t,
amíg csak az integrál értéke előírt mértékben pontosítható.
Bocs a hosszuert. Sikítsatok ha valami gázos.
--
Józsi
|
