5 g, Informatyka, Informatyka, Informatyka. Metody numeryczne, Kosma Z - Metody i algorytmy numeryczne [2009], Kosma Z

15 1.500000E-01 5.020501E-01 -7.229887E-06

16 1.600000E-01 5.168618E-01 -6.675536E-06

17 1.700000E-01 5.279421E-01 1.326847E-05

18 1.800000E-01 5.355462E-01 2.217025E-05

19 1.900000E-01 5.399292E-01 1.354511E-05

20 2.000000E-01 5.413460E-01 -4.885244E-06

21 2.100000E-01 5.400530E-01 -2.014654E-05

22 2.200000E-01 5.363109E-01 -2.363759E-05

23 2.300000E-01 5.303818E-01 -1.249445E-05

24 2.400000E-01 5.225275E-01 7.882037E-06

25 2.500000E-01 5.130102E-01 2.109110E-05

26 2.600000E-01 5.020790E-01 1.035045E-05

27 2.700000E-01 4.899333E-01 -5.095727E-06

28 2.800000E-01 4.767594E-01 -8.509673E-06

29 2.900000E-01 4.627439E-01 1.376860E-06

30 3.000000E-01 4.480733E-01 1.030785E-05

31 3.100000E-01 4.329220E-01 8.328657E-07

32 3.200000E-01 4.174159E-01 -1.091432E-05

33 3.300000E-01 4.016687E-01 -1.102860E-05

34 3.400000E-01 3.857943E-01 7.115909E-07

35 3.500000E-01 3.699064E-01 1.158983E-05

36 3.600000E-01 3.541091E-01 6.355500E-06

37 3.700000E-01 3.384682E-01 -3.127235E-06

38 3.800000E-01 3.230398E-01 -5.878178E-06

39 3.900000E-01 3.078801E-01 -5.869283E-07

40 4.000000E-01 2.930450E-01 5.193792E-06

................................................

Odchylenie kwadratowe: 5.661560E-07

5.7. Aproksymacja średniokwadratowa funkcji
dwóch zmiennych

Podobnie jak przy interpolacji funkcji wielu zmiennych najbardziej użyteczne praktycznie są metody aproksymowania funkcji wielu zmiennych oparte na zastosowaniu wielomianowych funkcji sklejanych.

Rozważymy aproksymację średniokwadratową funkcji dwóch zmiennych niezależnych określonej na dyskretnym zbiorze argumentów w obszarze prostokątnym
dwusześcienną funkcją sklejaną
przedstawioną przez znormalizowane B-funkcje sklejane (4.118). Uzmienniając współczynniki we wzorze (4.120)

otrzymamy

(5.110)

Dla funkcji zadanej na dyskretnym zbiorze argumentów minimalizujemy uogólniony funkcjonał (5.64)

(5.111)

Z warunków na minimum tego funkcjonału

wynika układ równań

(5.112)

gdzie:

Układ równań (5.112) jest rozwiązywany w programie 5.5, napisanym w języku Fortran ze względu na konieczność użycia tablic o dużych rozmiarach.

c {program 5.5}

program aproks

integer rozm

real*8 Bf,dx1,dx2,dy1,dy2,fun,ra,sm,xmin,xmax,ymin,ymax,x,y,z

real*8 S(-1:9,-1:9,-1:9,-1:9),mac(121,122),a(-1:9,-1:9),

* V(-1:9,-1:9),zw(0:30,0:30)

common/obsz/mac

c open (unit=1,file='pr_5_5.dan')

open (unit=4,file='pr_5_5.wyn')

write(*,*)'PROGRAM 5.5'

write(*,*)'Aproksymacja funkcji dwoch zmiennych.'

write(*,*)'Bikubiczna funkcja sklejana.'

write(*,*)'Rozmiary siatki wezlow - N, M: '

read(*,*)N,M

write(*,*)'Liczby linii kraty - K, L: '

read(*,*)K,L

write(*,*)'xmin, xmax, ymin, ymax: '

read(*,*)xmin,xmax,ymin,ymax

dx1=(xmax-xmin)/K

dy1=(ymax-ymin)/L

dx2=(xmax-xmin)/N

dy2=(ymax-ymin)/M

1010 format(2x,i3,4x,i3)

1020 format(2x,1pe13.6,2x,e13.6,2x,e13.6,2x,e13.6)

write(4,1010)K,L

write(4,1020)xmin,xmax,ymin,ymax

do 60 ii=-1,N+1

do 50 jj=-1,M+1

do 40 mi=-1,N+1

do 30 ni=-1,M+1

sm=0

do 20 i=0,K

x=xmin+i*dx1

do 10 j=0,L

y=ymin+j*dy1

sm=sm+Bf(ii,x,xmin,dx2)*Bf(jj,y,ymin,dy2)*

* Bf(mi,x,xmin,dx2)*Bf(ni,y,ymin,dy2)

10 continue

20 continue

S(ii,jj,mi,ni)=sm

30 continue

40 continue

50 continue

60 continue

do 100 mi=-1,N+1

do 90 ni=-1,M+1

sm=0

do 80 i=0,K

x=xmin+i*dx1

do 70 j=0,L

y=ymin+j*dy1

call random_number(ra)

c z=fun(x,y)+0.1*(0.5-ra)

z=fun(x,y)*(1+0.3*(0.5-ra))

c z=fun(x,y)

sm=sm+z*Bf(mi,x,xmin,dx2)*Bf(ni,y,ymin,dy2)

70 continue

80 continue

V(mi,ni)=sm

90 continue

100 continue

rozm=(N+3)*(M+3)

kk=0

do 120 mi=-1,N+1

do 120 ni=-1,M+1

ll=0

kk=kk+1

mac(kk,rozm+1)=V(mi,ni)

do 110 i=-1,N+1

do 110 j=-1,M+1

ll=ll+1

110 mac(kk,ll)=S(i,j,mi,ni)

120 continue

call elgaussa(rozm,1)

kk=0

do 130 mi=-1,N+1

do 130 ni=-1,M+1

kk=kk+1

130 a(mi,ni)=mac(kk,rozm+1)

odch=0

do 160 kk=0,K

x=xmin+kk*dx1

do 150 ll=0,L

y=ymin+ll*dy1

z=0

do 140 i=-1,N+1

do 140 j=-1,M+1

140 z=z+a(i,j)*Bf(i,x,xmin,dx2)*Bf(j,y,ymin,dy2)

zw(kk,ll)=z

odch=odch+(z-fun(x,y))*(z-fun(x,y))

150 continue

160 continue

write(4,1020)((zw(kk,ll),ll=0,L),kk=0,K)

write(*,*)'Odchylenie kwadratowe = ',odch

close (unit=1)

close (unit=4)

stop

end

real*8 function fun(x,y)

real*8 x,y,pi

pi=4*datan(1d0)

fun=dsin(2*dsqrt(x*x+y*y))

return

end

real*8 function Bf(i,x,a,h)

real*8 x,a,h,s,xx

s=0

xx=a+i*h

xx=dabs((x-xx)/h)

if(xx.ge.0.and.xx.le.1)s=3*xx*xx*xx-6*xx*xx+4

if(xx.gt.1.and.xx.le.2)then

s=2-xx

s=s*s*s

endif

Bf=s/6

return

end

subroutine elgaussa(n,r)

integer r,wm,w,k,p

real*8 a,d,m

common/obsz/a(121,122)

do 50 p=1,n-1

d=0d0

wm=p

do 10 w=p,n

m=dabs(a(w,p))

if(m.gt.d)then

d=m

wm=w

endif

10 continue

if(d.ne.0)then

do 20 k=p,n+r

m=a(p,k)

a(p,k)=-a(wm,k)

a(wm,k)=m

20 continue

do 40 w=p+1,n

m=a(w,p)/a(p,p)

do 30 k=n+r,p,-1

if(m.ne.0)a(w,k)=a(w,k)-a(p,k)*m

30 continue

40 continue

endif

50 continue

c d=1.0d0

c do 60 w=1,n

c 60 d=d*a(w,w)

c det=d

if(d.ne.0.and.r.gt.0)then

do 90 k=n+1,n+r

a(n,k)=a(n,k)/a(n,n)

do 80 w=n-1,1,-1

m=0d0

do 70 p=n,w+1,-1

70 m=m+a(w,p)*a(p,k)

a(w,k)=(a(w,k)-m)/a(w,w)

80 continue

90 continue

endif

return

end

Przy wykorzystaniu programu 5.5 aproksymowano funkcje (4.148) ÷ (4.151) określone na dyskretnym zbiorze argumentów, rozmieszczonych na płaszczyźnie
w węzłach prostokątnej siatki o rozmiarach
Jako węzły aproksymującej funkcji sklejanej (5.110) przyjęto węzły równomiernej prostokątnej siatki
dla Po wyznaczeniu współczynników funkcji (5.100) obliczono jej wartości na liniach wyjściowej kraty
Otrzymane wyniki były wczytywane przez program zamieszczony w podręczniku [15] i w ten sposób uzyskano obrazy aproksymowanych funkcji, pokazane na rysunkach 5.32 ÷ 5.35. Oprócz tego aproksymowano dyskretne zbiory wartości funkcji (4.148) ÷ (4.151) zaburzone funkcją losową o stałej amplitudzie - wynoszącej 0.05 (rys. 5.36 ÷ 5.39) oraz dyskretne zbiory wartości funkcji, będące iloczynami wartości funkcji zadanych (4.148) i (4.149) i funkcji losowej o amplitudzie 0.15 (rys. 5.40 - 5.41).