Аппроксимация таблично заданных функций

Министерство науки и  образования РФ

Новосибирский государственный  технический университет

 

Кафедра аэрогидродинамики

 

 

 

 

 

Отчет

по лабораторной работе №4

 

АППРОКСИМАЦИЯ ТАБЛИЧНО ЗАДАННЫХ ФУНКЦИЙ

 

 

 

 

 

Выполнил:

Студент гр. СА-01

Будко А.Н.

 

Проверил:

Гостеев Ю.А.

 

 

 

 

 

 

 

 

 

 

 

 

 

Новосибирск

2012

Задание:

В задаче аппроксимации (сглаживания) данных строится многочлен,

например, 3-й степени,

F(x) = a₁ + a₂x + a₃x² + a₃x³ , (1)

имеющий в узлах таблицы x_i минимальное отклонение от заданных

значений f_i. В i-й точке функция F(x) отклоняется от значения f_i на

величину f_i − F(x_i ) . Суммируя квадраты отклонений по всем точкам

i =1,2,..., N , находим суммарную невязку:

Потребуем, чтобы Φ(a₁,a₂,a₃, a₄)→ min. Это будет выполняться, если

Собирая коэффициенты при a_i , получим СЛАУ относительно искомых

коэффициентов:

или

где B – симметричная положительно определенная матрица. Решение (2)

можно получить с помощью  точного или приближенного методов  решения

СЛАУ (например, метода Гаусса или итерационного метода).

 

 

Текст программы:

program MNK

use Graflib

parameter(N=11)

dimension x(N),f(N),z(101),P(101),B(4,4),A(4),C(4)

open(1,file='tab.dat')

do i=1,N

read(1,*),x(i),f(i)

print*,x(i),f(i)

end do

B(1,:)=(/ 1.*N ,sum(x), sum(x**2),sum(x**3) /)

B(2,:)=(/ sum(x) ,sum(x**2),sum(x**3),sum(x**4) /)

B(3,:)=(/ sum(x**2),sum(x**3),sum(x**4),sum(x**5) /)

B(4,:)=(/ sum(x**3),sum(x**4),sum(x**5),sum(x**6) /)

print*,'B='

print '(4f10.3)',B

C =(/ sum(f),sum(f*x),sum(f*x**2),sum(f*x**3) /)

print*,'C=';print '(4f10.3)',C

call GAUSS (4,B,C,A)

print*,' A=';print*,A

do i=1,101;z(i)=x(1)+(i-1)*(x(N)-x(1))/100

end do

P=a(1)+a(2)*z+a(3)*z**2+a(4)*z**3

xgmin=minval(x);xgmax=maxval(x)

ygmin=minval(f)-0.5;ygmax=maxval(f)

call GrafInit(1,'MNK sample')

call Axis('X','Y',2)

call Plot(dble(x),dble(f),5,0,2)

call Plot(dble(z),dble(P),1,0,1)

end

 

Текст графического модуля:

module GrafLib

use dflib

real*8 xgmin,xgmax,ygmin,ygmax,zgmin,zgmax

common /gr1/ xgmin,xgmax,ygmin,ygmax,zgmin,zgmax

type (wxycoord) wxy

type(qwinfo)qww

contains

subroutine GrafInit(win,tit) ! инициализация графики

integer win; character(*) tit

qww.type=2;ii=setwsizeqq(QWIN$FRAMEWINDOW,qww)

open(win,file='user',title=tit)

call setviewport(10,10,810,810)

dxg=dabs(xgmax-xgmin); dyg=dabs(ygmax-ygmin)

xgmin1=xgmin-0.1*dxg; xgmax1=xgmax+0.1*dxg

ygmin1=ygmin-0.1*dyg; ygmax1=ygmax+0.1*dyg

ii=setwindow(.true.,xgmin1,ygmax1,xgmax1,ygmin1)

ncol=setcolor(15) ! задание текущего цвета (0...15)

ii=rectangle_w(3,xgmin1,ygmax1,xgmax1,ygmin1) ! заливка окна

if(initializefonts()<0) stop 'fonts not found!' ! инициализация шрифтов

if(setfont('t''Times New Roman Cyr''h30''b')==0) stop 'font not set!'

end subroutine GrafInit

subroutine Axis(xax,yax,icol) ! построение 2D-осей

character(*) xax,yax

ncol=setcolor(icol)

call moveto_w(xgmin,0.0_8,wxy); iz=lineto_w(xgmax,0.0_8)

call moveto_w(0.0_8,ygmin,wxy); iz=lineto_w(0.0_8,ygmax)

call moveto_w(xgmax,0.0_8,wxy); call outgtext(xax)

call moveto_w(0.0_8,ygmax,wxy); call outgtext(yax)

end subroutine Axis

subroutine Plot(x,y,icol,isl,im) ! построение 2D-кривой

real*8 x(:),y(:)

imin=lbound(x,dim=1); imax=ubound(x,dim=1)

ncol=setcolor(icol)

select case(im);case(1);mask=#FFFF;case(2);mask=#FF00;end select

call setlinestyle(mask)

call moveto_w(x(imin),y(imin),wxy)

do i=imin,imax

iz=lineto_w(x(i),y(i)); call sleepqq(isl)

end do

end subroutine Plot

end module GrafLib

 

Текст модуля Гаусса:

SUBROUTINE GAUSS (N,A1,b1,X)

! Метод Гаусса  для СЛАУ A1*X=b1

DIMENSION A1(N,N),A(N,N+1),b1(N),X(N),NI(N)

! формируем расширенную матрицу

do j=1,N;A(:,j)=A1(:,j);end do;A(:,N+1)=b1

DO J=1,N

NI(J)=0

END DO

NJ=N+1

DO IJ=1,N

B=0.

DO I=1,N

IF (NI(I)==0) THEN

DO J=1,N

C=ABS(A(I,J))

D=ABS(B)

IF (C>D) THEN

B=A(I,J)

K=I

L=J

END IF

ENDDO

END IF

END DO

NI(K)=L

A(K,L)=0.

D=ABS(B)

IF (D==0.0) THEN

PRINT *,'MATRIX IS SINGULAR!'

STOP

END IF

DO I=1,NJ

C=ABS(A(K,I))

IF (C/=0.0) THEN

S=A(K,I)/B

A(K,I)=S

DO J=1,N

A(J,I)=A(J,I)-A(J,L)*S

END DO

END IF

END DO

DO J=1,N

A(J,L)=0.

ENDDO

END DO

DO J=1,N

K=NI(J)

X(K)=A(J,N+1)

END DO

RETURN

END

function RESID(N,A,x,b)

! Невязка  решения СЛАУ R=max(abs(A*xb))

dimension A(N,N),x(N),b(N)

RESID=maxval(abs(matmul(A,x)-b))

end 

Результат программы: