| |
|
|
p.specht
| One of US-Volkswirt Hiroshi Akima 1978 vorgestelltes take action uses to the usual "linearen Interpolation" between Stützstellen Polynome third Grades. In its Berechnung from whom too unregelmäßig distributed - Stützstellen (Messwertetabelle) weg too The each benachbarten Stützstellen one as well as additional one gewogener rose one.
Window Title mkstr$(" ",26)+upper$("AKIMA-Interpolation through Kubischer Splines [H. Akima, 1978]")
' Portiert to XProfan11 (CL) CopyLeft 2014-08 by P. woodpecker, Vienna (Austria).
'{ No however geartete warranty. Use on alleinige menace the Anwenders!
'
' Übersetzungsbasis: Turbo Pascal manuals 1974; for Detailklärungen watts uses:
' F90/95: https://www.tacc.utexas.edu/documents/13601/162125/fortran_class.pdf
' F95: https://www.mrao.cam.ac.uk/~rachael/compphys/SelfStudyF95.pdf
'
' Pascal-DemoSource: https://jean-pierre.moreau.pagesperso-orange.fr/Pascal/akima_pas.txt
' Orig. 2D-Algorithmus: Hiroshi Akima,u.s. Department of Commerce, NTIA/ITS, F90-Version of 1995/08,
' published in "Hiroshi Akima, 'A Method of Bivariate Interpolation and Smooth Surface Fitting
' for Irregularly Distributed Data Points", ACM Transactions on Math. software Vol.4 No.2, June 1978
' pp. 148-159. copyright 1978, Association for Computing Machinery, Inc.
' Fortran-90: TOMS 22.yr,No3 (Sept. 1996) p.357ff, https://calgo.acm.org/760.gz (Abfrg. 2014-08)
' Test Data: https://cran.r-project.org/web/packages/akima/akima.pdf
'
'********************************************************
'* Program to demonstrate the Akima spline fitting *
'* of Function SIN(X) in double precision *
'* ---------------------------------------------------- *
'* Reference: BASIC Scientific Subroutines, Vol. II *
'* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. *
'* *
'* Pascal Version by J.-P. Moreau, Paris. *
'* (www.jpmoreau.fr) *
'* ---------------------------------------------------- *
'* SAMPLE RUN: *
'* *
'* Akima spline fitting of SIN(X): *
'* *
'* X SIN(X) HANDBOOK AKIMA INTERPOLATION ERROR *
'* ---------------------------------------------------- *
'* 0.00 0.0000000 0.0000000 0.0000000 *
'* 0.05 0.0499792 0.0500402 -0.0000610 *
'* 0.10 0.0998334 0.0998435 -0.0000101 *
'* 0.15 0.1494381 0.1494310 0.0000072 *
'* 0.20 0.1986693 0.1986459 0.0000235 *
'* 0.25 0.2474040 0.2474157 -0.0000118 *
'* 0.30 0.2955202 0.2955218 -0.0000016 *
'* 0.35 0.3428978 0.3428916 0.0000062 *
'* 0.40 0.3894183 0.3894265 -0.0000081 *
'* 0.45 0.4349655 0.4349655 0.0000000 *
'* 0.50 0.4794255 0.4794204 0.0000051 *
'* 0.55 0.5226872 0.5226894 -0.0000021 *
'* 0.60 0.5646425 0.5646493 -0.0000068 *
'* 0.65 0.6051864 0.6051821 0.0000043 *
'* 0.70 0.6442177 0.6442141 0.0000035 *
'* 0.75 0.6816388 0.6816405 -0.0000017 *
'* ---------------------------------------------------- *
'* *
'*******************************************************}
'}
PROC Interpol_Akima
'********************************************************
'* Akima Spline Fitting Subroutine *
'* -----------------------------------------------------*
'* The input table is X(i), Y(i), where Y(i) is the *
'* dependant variable. The interpolation point is xx, *
'* which is assumed to be in the interval of the table *
'* with at least one table value to the left, and three *
'* to the right. The interpolated returned value is yy. *
'* n is returned as on error checked (n=0 implies error). *
'* It is means assumed that the X(i) are in ascending *
'* order! *
'********************************************************
parameters i&
n&=1
' {special case xx=0}
if xx!=0: yy!=0:return:endif
' {Check to see if interpolation point is correct}
if (xx!<X![1]) or (xx!>=X![iv&-3]): n&=0:return:endif
X![0]=2*X![1]-X![2]
' {Calculate Akima coefficients, a and b}
whileloop iv&-1:i&=&Loop
' {Invoice values i to i+2}
XM![i&+2]=(Y![i&+1]-Y![i&])/(X![i&+1]-X![i&])
endwhile
XM![iv&+2]=2*XM![iv&+1]-XM![iv&]
XM![iv&+3]=2*XM![iv&+2]-XM![iv&+1]
XM![2]=2*XM![3]-XM![4]
XM![1]=2*XM![2]-XM![3]
whileloop iv&:i&=&Loop
a!=ABS(XM![i&+3]-XM![i&+2])
b!=ABS(XM![i&+1]-XM![i&])
if (a!+b!)<>0
Z![i&]=(a!*XM![i&+1]+b!*XM![i&+2])/(a!+b!)
else
Z![i&]=(XM![i&+2]+XM![i&+1])/2
endif
endwhile
' {Find relevant table interval}
i&=0
repeat
inc i&
until xx!<=X![i&]
dec i&
' {Begin interpolation}
b!=X![i&+1]-X![i&]
a!=xx!-X![i&]
yy!=Y![i&]+Z![i&]*a!+(3*XM![i&+2]-2*Z![i&]-Z![i&+1])*a!*a!/b!
yy!=yy!+(Z![i&]+Z![i&+1]-2*XM![i&+2])*a!*a!*a!/(b!*b!)
ENDPROC
'{ Main Program AKIMA-SINTEST
Windowstyle 24:font 0:set("decimals",8)
CLS'Window 0,0-%maxx,%maxy
var fmt$=" 0.0000000;-0.0000000"
var SIZE& = 14
declare X![SIZE&+1],Y![SIZE&+1],XM![SIZE&+4],Z![SIZE&+1]
declare a!,b!,xx!,yy!, i&,iv&,n&
iv&=14'{Number of points in table}
CLS
Font 2:Print "\n";tab(10);" Hiroshi Akima's Interpolation with kubischen Splines"
Font 0:print tab(10);" "+mkstr$("-",53)
print tab(10);" Base are 14 unregelmäßig distributed Stützwerte of/ one "
print tab(10);" Eigenbau-Sinusfunktion. Errechnet becomes the absolute "
Print tab(10);" Error (as difference to internen Sin(x)-function). "
print tab(10);" in the Programmtext are The Vergleichswerte to find! \n"
' -----------------------------------------------------------------
' Input y=Sine(x) table
' -----------------------------------------------------------------
' Sine table values from Handbook of mathematical functions
' by M. Abramowitz and I.A. Stegun, NBS, june 1964
' -----------------------------------------------------------------
X![1]=0.000:Y![1]=0.00000000 : X![2]=0.125:Y![2]=0.12467473
X![3]=0.217:Y![3]=0.21530095 : X![4]=0.299:Y![4]=0.29456472
X![5]=0.376:Y![5]=0.36720285 : X![6]=0.450:Y![6]=0.43496553
X![7]=0.520:Y![7]=0.49688014 : X![8]=0.589:Y![8]=0.55552980
X![9]=0.656:Y![9]=0.60995199 : X![10]=0.721:Y![10]=0.66013615
X![11]=0.7853981634 : Y![11]=0.7071067812
X![12]=0.849:Y![12]=0.75062005: X![13]=0.911:Y![13]=0.79011709
X![14]=0.972:Y![14]=0.82601466
' -----------------------------------------------------------------
print "\n";mkstr$("-",77)
font 2
print tab(3+3);"X";tab(%pos+15);"SIN(X)";\
tab(%pos+3);"AKIMA-Stützwertinterpolation";tab(%pos+2);"ABS.FEHLER"
font 0
print mkstr$("=",77)
' {main loop}
xx!=0
whileloop 16:i&=&Loop
Interpol_Akima(i&)
print tab(4);xx!;tab(%pos+6);SIN(xx!);tab(%pos+9);yy!;tab(%pos+11);stature$(fmt$,SIN(xx!) - yy!)
casenot i& mod 4:print
xx! = xx! + 0.05
endwhile
' {print footer}
print mkstr$("-",77)
beep
waitinput
END
'} of file akima.prf
|
|
|
| XProfan 11Computer: Gerät, daß es in Mikrosekunden erlaubt, 50.000 Fehler zu machen, zB 'daß' statt 'das'... | 05/11/21 ▲ |
|
|
|