INTERPOLACE KUBICKÝMI SPLINY
PROBLÉM
Jsou dány hodnoty funkce v bodech A.1<A.2<...<A.N. Odhadněte hodnotu funkce v bodě V z intervalu [A.1,A.N].
DEFINICE
Funkce S se nazývá kubický spline na intervalu [A.1,A.N] jestliže
- S je na [A.1,A.N] definována;
- S, první derivace S, druhá derivace S jsou spojité funkce na [A.1,A.N];
- pro body A. funkce S platí
A.1<A.2<...<A.N a S je polynom stupně <=3 na každém intervalu [A.J,A.Jp1], kde Jp1=J+1.
ALGORITMUS
Je dáno pole A.1,...,A.N a B.1,...,B.N. Pole B. obsahuje hodnoty funkce v bodech A., platí A.1<A.2<...<A.N. Jsou dány hodnoty první derivace funkce v bodech A.1 a A.N, které označíme B1A1 a B1AN. Podprogram SPLINE vypočte hodnoty pole B2. - druhé derivace interpolované funkce v zadaných bodech A.; funkce SPLINE je volána jen jednou. Hodnotu funkce v bodě V určí interpolační funkce SPLINT.
IMPLEMENTACE SPLINE
Jednotka: vnitřní podprogram
Globální proměnné: vstupní rostoucí posloupnost A.1,...,A.N, vstupní pole B. a výstupní pole B2., které obsahuje vypočítané hodnoty druhé derivace interpolované funkce v zadaných bodech.
Parametry: přirozené číslo N, hodnoty B1A1 a B1AN prvních derivací interpolované funkce v bodech A.1 and A.N
SPLINE: procedure expose A. B. B2.
N = ARG(1); B1A1 = ARG(2); B1AN = ARG(3)
if B1A1 > 0.99E+99 then do; B2.1=0; U.1=0; end
else do
B2.1 = -0.5
U.1=(3 / (A.2 - A.1)) *,
((B.2 - B.1) / (A.2 - A.1) - B1A1)
end
do J = 2 to N - 1
Jm1 = J - 1; Jp1 = J + 1
Sig = (A.J - A.Jm1) / (A.Jp1 - A.Jm1)
P = Sig * B2.Jm1 + 2
B2.J = (Sig - 1) / P
U.J = (B.Jp1 - B.J) / (A.Jp1 - A.J) -,
(B.J - B.Jm1) / (A.J - A.Jm1)
U.J = (6 * U.J / (A.Jp1 - A.Jm1) - Sig * U.Jm1) / P
end
Nm1 = N - 1
if B1AN > 0.99E+99 then do; Qn = 0; Un = 0; end
else do
Qn = 0.5;
Un = (3 / (A.N - A.Nm1)) *,
(B1AN - (B.N - B.Nm1) / (A.N - A.Nm1))
end
B2.N = (Un - Qn * U.Nm1) / (Qn * B2.Nm1 + 1)
do J = N - 1 to 1 by -1
Jp1 = J + 1
B2.J = B2.J * B2.Jp1 + U.J
end
return
|
IMPLEMENTACE SPLINT
Jednotka: vnitřní funkce
Globální proměnné: vstupní pole: rostoucí posloupnost A., pole B. a pole B2. vypočítané podprogramem SPLINE
Parametry: přirozené číslo N, hodnota V
Vrací: odhad hodnoty B.V
SPLINT: procedure expose A. B. B2.
N = ARG(1); V = ARG(2)
L = 1; H = N
do while (H - L) > 1
K = (H + L) % 2
if A.K > V then H = K
else L = K
end
Q = A.H - A.L
if Q = 0 then
call ERROR "SPLINT: Vstupní posloupnost musí být rostoucí"
X = (A.H - V) / Q; Y = (V - A.L) / Q
BV = (X**3 - X) * B2.L + (Y**3 - Y) * B2.H
BV = X * B.L + Y * B.H + BV * (Q**2) / 6
return BV
ERROR: say ARG(1); exit
|
PŘÍKLAD Následující program zobrazí 2.52386364
2.71270431
do J = 1 to 6; A.J = J; end
B.1 = 1.1; B.2 = 2.5; B.3 = 2.6
B.4 = 3.0; B.5 = 5.0; B.6 = 4.0
call SPLINE 6, 0, 0
say SPLINT(6, 3.5); say SPLINT(6, 3.8)
exit
SPLINE: procedure expose A. B. B2.
...
SPLINT: procedure expose A. B. B2.
...
|
SOUVISLOSTI
SPOLUAUTOR
Literatura
Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical Recipes in C : the art of scientific computing - 2nd ed. University Press, Cambridge, 1992
|
