Кубическая сплайн-интерполяция: как рассчитать второй сплайн S2 - PullRequest
1 голос
/ 22 октября 2011

У меня в Паскале есть задача вычислить кубический сплайн. Я сделал большую часть программы, и у меня все еще есть вопрос. Моя задача для функции y = f (x) найти сплайны S1 и S2, также найти
delta1 = max | f (x) -S (x) |, delta2 = max | f '(x) -S1' (x) |, delta3 = max | f (x) -S2 (x) |, delta4 = max | е '(х) -S2' (х) |

Какая формула здесь полезна для S2?

Вот мой код:

program spline_1_2;
uses crt;
const
k=1000;
var
input_file:text;
out1_file:text;
out2_file:text;
n,i:integer;
a_int,b_int,h:real;
delta_1, delta_2:real;
x, alpha, beta, y, a_coef, b_coef, c_coef, d_coef:array[0..k] of real;
S, S1, S2:array[0..k] of real;
A, B, C, F, z:array[0..k] of real;

Function f_given(x:real):real;
 begin
  f_given:=Cos((x*x)/3)-Exp(Sin(2*x));
 end;
Function f_prime(x:real):real;
 begin
  f_prime:=-Cos((x*x)/3)*(1/3)*2*x-Exp(Sin(2*x))*Cos(2*x)*2;
 end;

begin
clrscr;

Assign(input_file,'in.txt');
Reset(input_file);

Assign(out1_file,'out1.txt');
Rewrite(out1_file);

Assign(out2_file,'out2.txt');
Rewrite(out2_file);

while not EOF(input_file) do 
 begin
  readln(input_file,a_int,b_int);
  readln(input_file,h);
  readln(input_file,n);
end;  
 for i:=0 to n do
  begin
    x[i]:=a_int+(i-0.5)*h;
    f_given(x[i]); 
    y[i]:= f_given(x[i]);
    f_prime(x[i]);

   A[i]:= h;
  C[i]:= 4*h;
  B[i]:= h;

 end;

for i:=1 to n-1 do
begin 
 a_coef[i]:=y[i];
 F[i]:=((y[i+1]-2*y[i]+y[i-1])*3)/(h*h);
 A[1]:=0;
  C[n-1]:=0;
end;
alpha[1]:=-C[1]/B[1];
beta[1]:=F[1]/B[1];

  for i:=1 to n-1 do
 begin
  z[i]:=(A[i]*alpha[i-1] + C[i]);
   alpha[i]:= -B[i]/z[i];
  beta[i]:= ( F[i]-A[i]*beta[i-1])/z[i];
end;

c_coef[n-1]:=(F[i] - A[i] * beta[n-1])/(C[i] + A[i] * alpha[n-1]);

for i:=n-2 downto 1 do
begin
c_coef[i]:=alpha[i]*c_coef[i+1] + beta[i];
c_coef[0]:=0;
c_coef[n]:=0;

 end;



for i:=0 to n-1 do
 begin 
  d_coef[i]:=(c_coef[i+1]-c_coef[i])/(3*h);
 end;
  b_coef[0]:=f_prime(x[0]);
  b_coef[n-1]:=f_prime(x[n-1]);

for i:=1 to n-2 do
 begin 
   b_coef[i]:=(h/3)*(c_coef[i+1]+2*c_coef[i]) + (y[i+1] - y[i])/h;
 end;

 //our spline
for i:=1 to n do
 begin
   S[i]:= a_coef[i] + b_coef[i]*(x[i-1]-x[i]) + c_coef[i]*sqr(x[i-1]-x[i]) +    d_coef[i]*sqr(x[i-1]-x[i])*(x[i-1]-x[i]);
   S1[i]:=b_coef[i]+2*c_coef[i]*(x[i-1]-x[i])+3*d_coef[i]*(x[i-1]-x[i])*(x[i-1]-x[i]);
 end;

for i:=1 to n-1 do
 begin 
   Writeln(out1_file, x[i],'  ',y[i]:3:6,'  ',S[i]:3:6);
   delta_1:=abs(y[i]-S[i]);
   Writeln(out2_file,x[i],'  ',f_prime(x[i]):3:6,'  ', S1[i]:3:6);
   delta_2:=abs(f_prime(x[i])-S1[i]);
  end;
   Writeln(out1_file,'"delta1" = ',delta_1);
   Writeln(out2_file,'"delta2" = ',delta_2);
   Readln;  
 Close(input_file);
 Close(out1_file);
 Close(out2_file);
end.

1 Ответ

0 голосов
/ 20 марта 2014

скорее всего уже устарел, но в любом случае сначала некоторые заметки

1. поскольку я понимаю, что у вас есть f (x) и вы хотите создать из него SPLINE

  • одна кубическая кривая или болееобъединенные сплайны?
  • в чем разница или значение S1, S2?
  • так что же такое f (x)?Я вижу f_given и f_prime, какой это?

2. после вычисления SPLINE необходимо вычислить отклонения

  • между SPLINE и исходным f (x)
  • и между их производными

3. каков ваш ввод (a, h, n)?

  • Я предполагаю, что это начало интервала x
  • h - шаг по оси x между точками
  • и n - количество точек для выборки

Как ее решить:

1.get SPLINE

  • просто получить контрольные точкииз f (x) для 4 баллов
  • и вычислите коэффициенты SPLINE (ваши a, b, c, d), решив уравнения SPLINE
  • после этого:

    t = < 0.0 , 1.0 >
    x = x0 + (x1-x0)*t
    y = a + b*t + c*t*t + d*t*t*t
    y'=     b   + c*t   + d*t*t   
    
  • , где x0, x1 - диапазон x

  • , а t - параметр SPLINE
  • , если у вас есть больше соединенных сегментов SPLINE, вычисляйте каждый сегментпо отдельности
  • не забудьте правильно соединить их, поэтому, если у вас есть точки p0, p1, p2, ...
  • , последовательность вызовов SPLINE будет

    p0,p0,p0,p1 // 3x p0 to ensure that SPLINE starts at point p0
    p0,p0,p1,p2
    p0,p1,p2,p3 // here continue normally 
    p1,p2,p3,p4
    p2,p3,p4,p5
    ...
    p5,p6,p7,p8
    p6,p7,p8,p9 // at the end the last point is also 3x
    p7,p8,p9,p9
    p8,p9,p9,p9
    

2. вычислить различияations)

  • просто сделайте для x = x0 .. x1
  • и вычислите разницу
  • сохраните максимальный результат
...