Текст подпрограммы и версий ip04r_p.zip , ip04e_p.zip |
Тексты тестовых примеров tip04r_p.zip , tip04e_p.zip |
Бикубическая интерполяция функции двух переменных, заданной вместе с частными производными первого порядка и смешанной производной второго порядка в вершинах элементарного сеточного прямоугольника.
Пусть X1L, X1U, X2L и X2U следующим образом определяют вершины элементарного сеточного прямоугольника:
(X1L,X2U) | | (X1U,X2U) ______________________________________ | 4 3 | | (X1,X2) | | 1 2 | ______________________________________ (X1L,X2L) | | (X1U,X2L)
Эти вершины пронумерованы на рисунке цифрами, начиная с нижней левой вершины против часовой стрелки. Пусть в каждой из вершин известны значения интерполируемой функции
Y = f(x1, x2) , обеих ее частных производных первого порядка: Y1 = ∂f / ∂x1 , Y2 = ∂f / ∂x2 и смешанной производной Y12 = ∂2 f / ∂x1 ∂x2 .
Подпрограмма IP04R при помощи бикубической интерполяции вычисляет в точке (X1,X2), заданной внутри исходного элементарного сеточного прямоугольника, значение функции
AY = f (X1, X2) , обеих ее частных производных первого порядка: AY1 = ∂ f (X1, X2) / ∂x1 , AY2 = ∂ f (X1, X2) / ∂x2 и смешанной производной AY12 = ∂2 f (X1, X2) / ∂x1 ∂x2.
Используются следующие формулы:
4 4 AY = ∑ ∑ ci j t i-1 u j-1 ; i=1 j=1 4 4 AY1 = ∑ ∑ (i - 1) ci j t i-2 u j-1 ; i=1 j=1 4 4 AY2 = ∑ ∑ (j - 1) ci j t i-1 u j-2 ; i=1 j=1 4 4 AY12 = ∑ ∑ (i - 1) (j - 1) ci j t i-2 u j-2 ; i=1 j=1
Здесь t = (X1 - X1L) / (X1U - X1L), U = (X2 - X2L) / (X2U - X2L). Коэффициенты ci j ( i, j = 1, 2, 3, 4 ) вычисляются в подпрограмме IP04R.
Н.С.Бахвалов. Численные методы. Изд - во "Наука", 1973.
procedure IP04R(var Y :Array of Real; var Y1 :Array of Real; var Y2 :Array of Real; var Y12 :Array of Real; X1L :Real; X1U :Real; X2L :Real; X2U :Real; X1 :Real; X2 :Real; var AY :Real; var AY1 :Real; var AY2 :Real; var AY12 :Real; var C :Array of Real);
Параметры
Y, Y1 - Y2, Y12 | вещественные векторы длины 4, содержащие соответственно значения интерполируемой функции, обеих ее частных производных первого порядка и смешанной производной в вершинах заданного элементарного сеточного прямоугольника, пронумерованных начиная с нижней левой вершины против часовой стрелки; |
X1L - X1U | левая и правая координаты вершин элементарного сеточного прямоугольника по первой переменной (тип: вещественный); |
X2L - X2U | нижняя и верхняя координаты вершин элементарного сеточного прямоугольника по второй переменной (тип: вещественный); |
X1, X2 - | координаты точки, заданной внутри исходного элементарного сеточного прямоугольника, в которой выполняется интерполяция (тип: вещественный); |
AY, AY1 - AY2, AY12 | вещественные переменные, значения которых полагаются равными соответственно значениям в точке (X1, X2) интерполируемой функции, ее частных производных первого порядка и смешанной производной; |
C - | вещественный двумерный массив размеров 4 на 4, в котором помещаются коэффициенты, используемые при бикубической интерполяции. |
Версии
IP04E - | бикубическая интерполяция функции двух переменных в режиме расширенной (Extended) точности; при этом все формальные параметры должны иметь тип Extended. |
Вызываемые подпрограммы: нет
Замечания по использованию: нет
Unit tip04r_p; interface uses SysUtils, Math, { Delphi } Lstruct, Lfunc, UtRes_p, IP04R_p; function tip04r: String; implementation function tip04r: String; var _i :Integer; X1L,X1U,X2L,X2U,X1,X2,TY,TY1,TY2,TY12,AY1,AY2,AY12,AY :Real; Y :Array [0..3] of Real; Y1 :Array [0..3] of Real; Y2 :Array [0..3] of Real; Y12 :Array [0..3] of Real; C :Array [0..15] of Real; begin Result := ''; { результат функции } X1L := 0.0; X1U := 0.1; X2L := 0.0; X2U := 0.15; Y[0] := Sin(X1L)*Sin(X2L); Y[1] := Sin(X1U)*Sin(X2L); Y[2] := Sin(X1U)*Sin(X2U); Y[3] := Sin(X1L)*Sin(X2U); Y1[0] := Cos(X1L)*Sin(X2L); Y1[1] := Cos(X1U)*Sin(X2L); Y1[2] := Cos(X1U)*Sin(X2U); Y1[3] := Cos(X1L)*Sin(X2U); Y2[0] := Sin(X1L)*Cos(X2L); Y2[1] := Sin(X1U)*Cos(X2L); Y2[2] := Sin(X1U)*Cos(X2U); Y2[3] := Sin(X1L)*Cos(X2U); Y12[0] := Cos(X1L)*Cos(X2L); Y12[1] := Cos(X1U)*Cos(X2L); Y12[2] := Cos(X1U)*Cos(X2U); Y12[3] := Cos(X1L)*Cos(X2U); X1 := 0.06; X2 := 0.07; IP04R(Y,Y1,Y2,Y12,X1L,X1U,X2L,X2U,X1,X2, AY,AY1,AY2,AY12,C); Result := Result + #$0D#$0A; for _i:=0 to 15 do begin Result := Result + Format('%20.16f ',[C[_i]]); if ( ((_i+1) mod 3)=0 ) then Result := Result + #$0D#$0A; end; Result := Result + #$0D#$0A; TY := Sin(X1)*Sin(X2); TY1 := Cos(X1)*Sin(X2); TY2 := Sin(X1)*Cos(X2); TY12 := Cos(X1)*Cos(X2); Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f %20.16f %20.16f %20.16f %20.16f ', [AY,AY1,AY2,AY12,TY,TY1,TY2,TY12]) + #$0D#$0A; UtRes('tip04r',Result); { вывод результатов в файл tip04r.res } exit; end; end. Результаты: AY = 0.4194o5E-02 AY1 = 0.698169E-01 AY2 = 0.598171E-01 AY12 = 0.995755