Kubinis splainas
1. Matematinis uždavinio formulavimas.
Duota funkcijos y = f(x) reikšmių lentelė (xi, yi) , i = 1, 2, . ,N. Rasti kubinį splainą y = S3i(x) , tenkinantį Lagranžo interpoliavimo sąlygą: S3i(x) = yi , i = 1, 2, . ,N.
x0=a x1 x2 . xn=b
y0 y1 y2 . yn
Reikia paskaičiuoti y-o reikšmę, pagal bet kokią x-o reikšmę iš intervalo [a, b], naudojant sudarytąjį kubinį splainą.
PVZ. Duota štai tokia lentelė.
xi 0 0.2 0.4 0.6 0.8 1 1.2
yI 1.2 4 0.8 2.5 2 3 1.5
Šio kubinio splaino grafikas:
Pagal lentelės duomenis sudarome kubinį ssplainą ir apskaičiuojame jo reikšmę taške x.
2. Teorinė dalis.
Splaino užrašymo formos.
Literatūroje nurodomos keturios splainų užrašymo formos, iš kurių tik dvi naudojamos praktiniuose skaičiavimuose, tai:
1) kiekviename intervale [xi, xi+1] nurodoma polinomo išraiška t.y., dalimis polinominė forma;
2) splainas užrašomas tiesine B – splainų kombinacija.
Splaino apibrėžimas.
Sakykime, kad intervale [a, b] yra apibrėžta funkcija y = f(x), ir tame intervale yra žinoma n+1 reikšmė.
x0=a x1 x2 . xn=b
y0 y1 y2 . yn
Funkcija S(x), kuri, pirma, kiekviename intervale [xi, xi+1] yra m – ojo llaipsnio daugianaris:
Si(x) = aixm + ai1xm-1 + . + aim-1x + aim , antra , tenkina interpoliavimo sąlygą S(xi)=f(xi), i = 0, 1, . ,n, ir, trečia, kurios išvestinės iki (m-1)-osios eilės imtinai yra tolydžios intervale [a, b], vadinama mm-osios eilės splainu.
Splainai gali būti tiesiniai, kvadratiniai, bet dažniausiai naudojami kubiniai splainai.
Kubiniai splainai.
Kubinis splainas S3(x) yra du kartus tolydžiai diferencijuojama funkcija – kubinis daugianaris kiekviename intervale xi x xi+1:
S3i(x) = aix3 + bix2 + cix + di , i = 0, 1, . ,n-1.
Iš viso kubinis splainas turi 4N koeficientų, taigi turime sudaryti 4N lygčių sistemą šiems koeficientams apskaičiuoti. Iš interpoliavimo sąlygos gausime 2N lygčių:
S3i(xi) = yi,
S3i(xi+1) = yi+1.
Iš išvestinių tolydumo sąlygos – dar (2N – 2) lygtis:
S’3i(xi) = S’3i-1(xi),
S”3i(xi) = S”3i-1(xi).
Taigi dar turime sudaryti dvi lygtis. Paprastai tos lygtys gaunamos iš kraštinių sąlygų.
Kraštinės sąlygos.
Praktiškai taikant kubinius splainus, reikia tinkamai parinkti kraštines sąlygas. Tai ypač svarbu, jei splainas turi gerai sutapti su aproksimuojamąja funkcija interpoliavimo intervalo ggalinėse srityse. Kubinių splainų kraštinės sąlygos didžiausią įtaką splaino formai turi interpoliavimo intervalo galinėse srityse ir turi nedidelę įtaką kitose interpoliavimo intervalo srityse. Tai reiškia, kad kraštines sąlygas galima parinkti, remiantis lokalinėmis aproksimuojamosios funkcijos savybėmis.
Jei žinomos f’(x0) ir f’(xN) arba f”(x0) ir f”(xN) reikšmės, tai tikslinga reikalauti, kad splaino S3i(x) atitinkamos eilės išvestinės taškuose x0 ir xN sutaptų su funkcijos f(x) tos pačios eilės išvestinėmis. Jei yra galimybė rinktis tarp pirmosios ir antrosios eilės išvestinių, pirmenybė teiktina pirmosios eilės iišvestinėms.
Tačiau dažniausiai yra keliamos sąlygos antrosios eilės išvestinėms:
S”30(x0) = 0,
S”3N(xN) = 0.
Šios sąlygos vadinamos natūraliosiomis, o splainas, tenkinantis natūraliąsias sąlygas – natūraliuoju. Fizinė šių sąlygų prasmė – splaino galai gali laisvai sukiotis apie ašį, esančią pradiniame (ar galiniame) mazge.
Pervertimas į triįstrižaininę lygčių sistemą.
Sudarytąją sistemą galima spręsti tiesiogiai, tačiau mes ją pertvarkysime taip, kad ji smarkiai supaprastėtų – paversime triįstrižaine lygčių sistema, kurioje yra (N+1) lygtis.
Kadangi S”3i(x) yra tiesinė funkcija, o S”3(x) yra tolydi, tai splaino antrąją išvestinę galime aproksimuoti tiesiniu splainu:
S”3i(x) = gi + (x – xi)(gi+1 – gi)/hi ,
čia gi ir gi+1 yra nauji nežinomieji,
gi = S”3i(xi).
hi = (xi – xi-1) – žingsnis tarp artimiausių mazgų.
Suintegruokime šią lygybę:
S’3i(x) = ei + gi(x – xi) + (gi+1 – gi)(x – xi)2/2hi ,
čia mes apibrėžėme dar vieną naują nežinomąjį S’3i(xi) = ei. Į šią lygtį įrašę x = xi+1, remdamiesi splaino išvestinės tolydumu, gausime:
ei+1 = ei + (gi+1 – gi)hi/2.
Dar karą suintegruokime tą pačią lygybę:
S3i(x) = yi + ei(x – xi) + gi(x – xi)2/2 + (gi+1 – gi)(x – xi)3/6hi .
Įrašykime į ją x = xi+1 ir išreikšime ei:
ei = (yi+1 – yi)/hi – gi+1hi/6 – gihi/3 .
nežinomojo ei+1 išraišką perrašykime i-ajame taške:
ei = ei-1 + (gi – gi-1)hi-1/2 == (yi – yi-1)/hi-1 – gihi-1/6 – gi-1hi-1/3 + (gi – gi-1)hi-1/2 .
Šių dviejų lygybių kairiosios pusės yra lygios, todėl ir dešiniosios pusės lygios:
(yi+1 – yi)/hi – gi+1hi/6 – gihi/3 = (yi – yi-1)/hi-1 – gihi-1/6 – gi-1hi-1/3 .
Pertvarkę lygtį, gausime:
6( (yi+1 – yi)/hi – (yi – yi-1)/hi-1) = gi-1hi-1 + 2gi(hi + hi-1) + gi+1hi .
Prijungę papildomąsias sąlygas, gausime triįstrižainę lygčių sistemą, kurioje yra (N–1) lygtis:
gi-1hi-1 + 2gi(hi + hi-1) + gi+1hi = 6( (yi+1 – yi)/hi – (yi – yi-1)/hi-1) , g0 = 0, gN = 0.
Šios sistemos sprendiniai yra nežinomieji gi. Juos radę, išreikštinėmis formulėmis apskaičiuosime koeficientus ei ir galėsime sudaryti splaino lygtį:
S3i(x) = yi + ei(x – xi) + Gi(x – xi)2 + Hi(x – xi)3 , i = 0, 1, . , N-1.
Čia pažymėjome Gi = gi/2, Hi = (gi+1 – gi)/6hi.
Triįstrižaininės lygčių sistemos sprendimas.
Triįstrižaininė lygčių sistema – kai lygčių sistemoje yra tik trys nenuliniai koeficientai: pagrindinėje įstrižainėje ir abiejose gretutinėse. Triįstrižaininės lygčių sistemos dažniausiai sprendžiamos modifikuotu Gauso metodu, vadinamu perkelties metodu.
Ji atrodo taip:
b1g1 + c1g2 = d1
a2g1 + b2g2 + c2g3 = d2
…
aigi-1 + bigi + cigi+1 = di
…
angn-1 + bngn = dn
Ją sprendžiant, įvesim pasižymėjimą: C1 = –c1/b1 , D1 = d1/b1 , Ci = -ci/(bi + Ci-1ai ) , Di = (di – Di-1ai)/(bi + Ci-1ai) , i = 2, 3, . ,n-1 .Paskutinėje, n-ojoje, lygtyje cn = 0, todėl Cn = 0 ir gn = Dn. Žinant gn reikšmę, kiti nežinomieji nuosekliai apskaičiuojami iš grįžtamosios formulės: gi = Ci gi+1 + Di .
Jei pagrindinė įstrižainė yra vyraujamoji, t. y. ai + ci≤bi, i = 1, 2, . ,n, ir bent su vienu i galioja griežta nelygybė, tai dalyba iš nulio ar labai mažo skaičiaus perkelties metodo eigoje negalima. Bet tada nesikaups apvalinimo paklaidos.
Taigi, perkelties algoritmas yra dviejų etapų: pirmiausia apskaičiuojami perkelties koeficientai Ci, Di, i = 2, 3, . ,n, o po to atbuline tvarka randamos nežinomųjų gi reikšmės.
Pirmuoju perkelties etapu skaičiuojant koeficientus Ci , Di , atliekama 4n-3 daugybų (su dalybomis) ir 2n-2 sudėčių (su atimtimis). Antruoju etapu skaičiuojant nežinimųjų gi reikšmes atliekama n-1 daugyba ir n-1 sudėtis. Taigi, sprendinys randamas atlikus 5n-4 daugybas ir 3n-3 sudėtis. Jei tarsime, kad kompiuteris daugina ir sudeda tuo pačiu greičiu, o n – didelis skaičius, tai perkelties metodo skaičiavimų apimtis proporcinga 8n. Tokia lygčių sistema išsprendžiama (1/12)n2 kartų greičiau negu Gauso metodu.
Splaino lygčių pavyzdys.
Sudarysime pagal lentelės
duomenis kubinį splainą:
xi 0 0.2 0.4 0.6 0.8 1 1.2
yi 1.2 4 0.8 2.5 2 3 1.5
Splaino lygtys atrodys taip:
-251.59×3 + 24.063x + 1.2 0 ≤ x ≤ 0.2 507.93(x – 0.2)3 – 150.95(x – 0.2)2 – 6.127(x – 0.2) + 4 0.2 ≤ x ≤ 0.4
417.64(x – 0.4)3 + 153.81(x – 0.4)2 – 5.556(x – 0.4) + 0.8 0.4 ≤ x ≤0.6
S3(x) 275.14(x – 0.6)3 – 96.78(x – 0.6)2 + 5.85(x – 0.6) + 2.5 0.6 ≤ x ≤0.8
-220.43(x – 0.8)3 + 68.31(x – 0.8)2 +0.156(x – 0.8) + 2 0.8 ≤ x ≤ 1
106.59(x – 1)3 + 63.95(x –1)2 + 1.027(x – 1) +3 1 ≤ x ≤ 1.2
Šio kubinio splaino grafikas pateiktas 1 skyriuje ‘Matematinio uždavinio formulavimas’.
Tikslumo įvertis.
Tikslumas įvertinamas dviem būdais. Jais kartu įvertinamas ir splaino pirmosios bei antrosios išvestinės tikslumas.
1) Teorema
Jei funkcija y(x)nyra keturis kartus tolydžiai difencijuojama intervale [a, b], t.y. y(x) C4[a, b], tai kubinio splaino S3(x) ir jo išvestinių S’3(x), S”3(x) paklaida įvertinama šia nelygybe:
y(k)(x) – S(k)3(x)≤ M4h2-k(h2 + max( y”(a) – S”3(a) ,, y”(b) – S”(b) )) , k = 0, 1, 2 .
Iš šios teoremos matyti, jei konstruojant kubinį splainą kraštinės sąlygos sutampa su tiksliomis interpoliuojamos funkcijos antrosios išvestinės y”(x) reikšmėmis, tai teoremos įvertis yra toks:
y(k)(x) – S(k)3(x)≤ M4h4-k .
Jei ssplaino kraštinės sąlygos nesutampa su y”(x) reikšmėmis, pavyzgžiui, pasirinkus natūraliąsias kraštines sąlygas, paklaidos įvertis yra toks:
y(k)(x) – S(k)3(x)≤ Mh2-k .
Konvergavimo pavyzdys.
Aproksimuosime funkciją y = ex kubiniais splainais intervale [0, 1] dviem atvejais: su tiksliosiomis ir natūraliosiomis kraštinėmis sąlygomis. Apskaičiuosime interpoliavimo paklaidą su skirtingu interpoliavimo mazgų skaičiumi N. 1 lentelėje pateikiama didžiausioji funkcijos ir išvestinių aproksimavimo paklaida su tiksliosiomis sąlygomis, 2 lentelėje – su natūraliosiomis kraštinėmis sąlygomis. Paklaidas pažymėjome:
║Z(k)║ = max y(k)(x) – S(k)3(x) , k = 0, 1, 2.
1 lentelė Kubinio splaino su tiksliosiomis kraštinėmis sąlygomis paklaida.
N 6 11 21 41
║Z║
║Z’║
║Z”║ 0.2675*10-4
0.4989*10-3
0.9817*10-2 0.1708*10-5
0.6386*10-4
0.2656*10-2 0.1079*10-6
0.8079*10-5
0.6904*10-3 0.6779*10-8
0.1016*10-5
0.1760*10-3
2 lentelė Kubinio splaino su natūraliosiomis kraštinėmis sąlygomis paklaida.
N 6 11 21 41
║Z║
║Z’║ 0.5257*10-2
0.1566 0.1317*10-2
0.0784 0.3294*10-3
0.0392 0.8239*10-4
0.0196
Matyti, kad du kartus padidinus mazgų skaičių, funkcijos aproksimavimo paklaida sumažėja apie 16 kartų, ppirmosios išvestinės – 8 kartus, o antrosios – 4 kartus. Tuo atveju, kai kraštinės sąlygos nėra aproksimuojamos tiksliai, tai paididnus du kartus mazgų skaičių funkcijos aproksimavimo paklaida sumažėja 4 kartus, o pirmosios išvestinės – 2 kartus.
2) Tarkime, kad segmente [a, b] yra duotas pastovaus žingsnio tinklelis, t. y.: a = x0 < x1 < . < xN = b ir
hi = xi+1 – xi = h , kai visi i = 1, . ,n. Apibrėšime funkcijos tolydumo modulį.
Tarkime, ssegmente [a, b] nusakyta funkcija (x). Tada dydis (h, ) = sup (x1) – (x2) vadinamas funkcijos (x) tolydumo moduliu, kuris parodo didžiausią pagal modulį funkcijos (x) pokytį ilgio h atkarpoje, kai h priklauso segmentui [a, b].
Tarkime, kad segmente [a, b] nusakyta tolydžioji funkcija f(x), turinti tolydžiasias pirmosios ir antrosios eilės išvestines. Šią funkciją interpoliuokime kubiniu splainu S3i(x), kurį nusako minėtas tinklelis ir kuris tenkina kraštinės sąlygas S’3i(x0) = f’(x0) , S’3i(xN) = f’(xN). Tada teisingi tokie rezultatai:
f”(x) – S3i(x) ≤ 5(h, f”),
f’(x) – S3i(x) ≤ 5h(h, f”),
f(x) – S3i(x) ≤ 5/2 h2 (h, f”).
Didėjant funkcijos glodumo laipsniui, t. y. kai funkcija f(x) turi aukštesnės nei 2 eilės tolydžiųjų išvestinių, konvergavimo greitis didėja. Tačiau negalima pasiekti didesnio konvergavimo greičio nei O(h4).
Kaip kubinis splainas aproksimuoja į funkciją cosx,kaip tai priklauso nuo mazgų skaičiaus parodysime grafiniu būdu.
Mėlyna linija pavaizduotas cosx grafikas.
Raudona linija pavaizduotas kubinio splaino grafikas su 5 mazgais.Šis kubinis spalinas sudaryta pagal lentelę:
Xi 0 /4 /2 3/4
Yi 1 0.707 0 -0.707 -1
Geltona linija pavaizduotas kubinio splaino grafikas su 9 mazgais.Šis kubinis splainas sudarytas pagal lentelę:
Xi 0 /8 /4 3/8 /2 5/8 3/4 7/8
Yi 1 0.924 0.707 0.383 0 -0.383 -0.707 -0.924 -1
Išvada: Iš viso to seka,kad didinant mazgų skaičių kubinis splainas greičiau aproksimuoja į duotąją funkciją,šiuo atveju į cosx funkciją.
3. Programos sudarymas.
Pagal lentelės duomenis sudarysime kubinį splainą ir apskaičiuosime jo reikšmę taške xx=0,4:
xI 0 0.2 0.4 0.6 0.8 1 1.2
yI 1.2 4 0.8 2.5 2 3 1.5
Programa Spline iš lentelės (x,y) sudaro kubinį splainą ir pagal jį apskaičiuoja funkcijos reikšmę norimame taške
program Spline
c Kintamuju aprasai
real xpr,sp
real x(10)
real y(10),sk(10),h(10),a(10),aa(10)
real b(10),c(10),g(10),dk(10),ck(10),e(10)
real gd(10),hd(10)
c Atidaromas failas i kuri bus rasomi rezultatai
Open(UNIT=9,FILE=’beavis.dat’)
c Duomenu ivedimas
10 write (*,*)’Iveskite x ir y reiksmiu skaiciu:’
read (*,*) n
if (n<=3) then ! Tikriname kad,ivestu reiksmiu butu reikiamas skaicius
write(*,*)’Tasku skaicius turi buti ne mazesnis uz 3 ‘
goto 10
end if
write(*,*)’Iveskite x reiksmes didejimo tvarka:’
do i=1,n
read (*,*) x(i)
end do
write(*,*)’Iveskite atitinkamas y reiksmes:’
read (*,*)(y(i),i=1,n)
c Funkcine dalis
c Zingsnio skaiciavimas
do i=1,n
h(i)=x(i)-x(i-1)
end do
do i=2,n-1
sk(i)=6*(((y(i+1)-y(i))/(h(i+1)))-((y(i)-y(i-1))/(h(i))))
end do
c Koeficiento a skaiciavimas
do i=1,n-1
a(i)=h(i)
end do
c Koeficiento b skaiciavimas
do i=1,n
b(i)=2*(h(i+1)+h(i))
end do
c Koeficiento c skaicivimas
do i=2,n
c(i)=h(i+1)
end do
c Ciklai triistrizaininei lygciu sistemai spresti ir g surasti
do i=2,n-1
ck(i)=(-c(i)/(b(i)+(ck(i-1)*a(i))))
dk(i)=(sk(i)-(dk(i-1)*a(i)))/(b(i)+(ck(i-1)*a(i)))
g(i)=dk(i)
end do
do i=n-1,2,-1
g(i)=ck(i)*g(i+1)+dk(i)
end do
c Ciklai kitu koeficientu radimui
write(9,*)’Koeficientai:’
write(9,*)’ yi ei Gi Hi’
do i=1,n-1
aa(i)=y(i)
e(i)=(y(i+1)-y(i))/h(i+1)-(g(i+1)*h(i+1))/6-(g(i)*h(i+1))/3
gd(i)=g(i)/2
hd(i)=(g(i+1)-g(i))/(h(i+1)*6)
write(9,*)’ ‘,aa(i),’ ‘,e(i),’ ‘,gd(i),’ ‘,hd(i) ! Koeficientu isvedimas i rezultatu faila
end do
c Ivedimas reiksmes x pagal kuria skaiciuosime kubini splaina
20 write(*,*)’Iveskite x reiksme:’
read (*,*) xpr
i=-1
c Ciklas apskaiciuoja kubini splaina su ivesta reiksme
do j=1,n-1
if((xpr.ge.x(j)).and.(xpr.le.x(j+1)))then ! Tikriname,kad ivesta x reiksme priklausytu anksciau sudarytam intervalui
i=j
sp=y(i)+(e(i)*(xpr-x(i)))+(gd(i)*(xpr-x(i))*(xpr-x(i))) ! Kubinio splaino apskaiciavimas
*+(hd(i)*(xpr-x(i))*(xpr-x(i))*(xpr-x(i)))
end if
end do
write(9,*)’x = ‘,xpr
write(9,*)’f(x)=’,sp
if (i.eq.-1) then ! Jei x nepatenka i intervala
write(*,*)’X turi buti vidinis taskas:’
goto 20
end if
end
4. Programos testavimas.
Koeficientai:
yi ei Gi Hi
1.200000 24.063460 0.000000E+00 -251.586500
4.000000 -6.126922 -150.951900 507.932600
8.000000E-01 -5.555770 153.807700 -417.644100
2.500000 5.849998 -96.778830 275.144200
2.000000 1.557699E-01 68.307690 -220.432700
3.000000 1.026925 -63.951910 106.586500
x = 4.000000E-01
f(x) = 8.000000E-01
5. Kontrolinio uždavinio sprendimas ir rezultatų pateikimas.
Kontorolinio uždavinio duomenys pateikti lentele.
Xi 0.030 0.085 0.261 0.270 0.451 0.577
Yi 1.020 1.057 1.172 1.172 1.290 1.364
Pagal lentelės duomenis sudarysime kubinį spaliną.Ir apskaičiuosime f(x) reikšmę mums duotame taške x=0.05.
Rezultatai:
Koeficientai:
yi ei Gi Hi
1.020000 6.765283E-01 0.000000E+00 -1.256187
1.057000 6.651284E-01 -2.072708E-01 7.993408E-01
1.172000 6.664503E-01 2.147811E-01 -21.301680
1.178000 6.651400E-01 -3.603645E-01 5.760039E-01
1.290000 5.912994E-01 -4.759436E-02 1.259110E-01
x = 5.000000E-02
f(x) = 1.033520
Kontrolinės užduoties kubinio splaino grafikas:
6. Išvados
