[octave] Полином Ньютона. Выход за пределы массива.

15.10.2011 23:35 #199882

Здравствуйте.

В лабораторной по вычислительной математике требуется интерполировать функцию. В том числе полиномом Ньютона.
Пишу:

text

C = [y];
for j=2:16
  tmp = [];
  for i=1:16
    if i > j
      ch = C(i)(j-1) - C(i-1)(j-1);
      zn = x(i) - x(i-j);
      c = ch / zn;
      tmp = [tmp; c];
    else
      tmp = [tmp; 0];
    end;
  end;
  C = [C, tmp];
end;

Где:

text

Y = [25, 27, 25, 26, 22, 23. 23. 20, 20, 21, 27, 24, 28, 27,
     25, 27, 28, 30, 32, 34, 32, 33, 33, 28, 29, 19, 18, 21,
     23, 22, 28];
X = 0:1:30;
y = [];
for i=1:2:31
  y = [y; Y(i)];
end;
x = 0:2:30;

На первый взгляд вроде все правильно, но octave мне выдает следующее:
[console]error: A(I): Index exceeds matrix dimension.[/console]

Что я делаю не так?

Заранее спасибо.

16.10.2011 00:27 #199889

Osanve написал:
C = [y];

надо объявить матрицу: C = [[y]];

Osanve написал:
ch = C(i)(j-1) - C(i-1)(j-1);

матрица индексируется так: ch = C(i, j-1) - C(i-1, j-1);

16.10.2011 00:28 #199891

Osanve написал:
C = [y];

надо объявить матрицу: C = [[y]];

Osanve написал:
ch = C(i)(j-1) - C(i-1)(j-1);

матрица индексируется так: ch = C(i, j-1) - C(i-1, j-1);

16.10.2011 01:23 #199896

Fat-Zer, спасибо, помогло.

Сейчас никак не могу понять, почему на графике интерполированную функцию выводит как прямую линию.

text

C = [[y]];
for j=2:16
  tmp = [];
  for i=1:16
    if i > j
      ch = C(i, j-1) - C(i-1, j-1);
      zn = x(i) - x(i-j);
      c = ch / zn;
      tmp = [tmp; c];
    else
      tmp = [tmp; 0];
    end;
  end;
  C = [C, tmp];
end;
F = [];
t = 0:0.1:30;
for k=0:0.1:30
  sum = 0;
  for i=1:16
    pr = 1;
    for j=1:16
      if i > j
        pr = pr * (k - x(j));
      end;
    end;
    sum = sum + (C(i, i) * pr);
  end;
  F = [F, sum];
end;
plot(X, Y, t, F)

В то же время при использовании канонического полинома - все нормально.

text

A = [];
for i=1:16
  tmp = [];
  for j=0:15
    tmp = [tmp, x(i)^j];
  end;
  A = [A; tmp];
end;
a = inv(A) * y;
t = 0:0.1:30;
P = [];
for i=0:0.1:30
  tmp = 0;
  for j=1:16
    tmp = tmp + (a(j) * (i^(j-1)));
  end;
  P = [P, tmp];
end;
plot(X, Y, t, P)

Заранее спасибо.