Построение интерполяционного сплайна степени 1

[Alexeybk5]

_Ivana в сообщении #638942 писал(а):

В чем могут быть сложности? Исходная функция на отрезке есть, задаете сетку, строите сплайн, анализируете модуль отклонения - хоть график его стройте, хоть максимум ищите, хоть анализируйте график этого максимума в зависимости от частоты сетки...

Вы правы, теоретически проблем нет, но как это в вольфраме реализовать?....
Кстати, ведь по сути, то что Вы сказали, я и пытался сделать. Функция sg1
должна на входе получаеть значения сетки разбиения и выводить для каждого x значения сплайна в этой точке.

Сложность в том, что функция не хочет печатать таблицу со значениями сплайна, в этих точках.
Т.е. в X[n] мы строим сетку, а в sg1 строим "подсетку" в точках которых находим значения сплайна sj, и эти значения выписываем в таблицу. После чего уже можем модулем разности в соответствующих подседках между значением сплайна и функции искать ошибку апроксимации...

Т.е. надо найти разность не на отрезке, а на совокупности отрезков. Т.е. Экспоненциальную функцию апроксимировать на одной прямой, я совокупностью прямых.

-- 02.11.2012, 00:50 --

ИСН в сообщении #638944 писал(а):

Перебором находить, между какими натуральными числами находится x/h - это, конечно, тот ещё прикол. Но ладно. Хуже другое: я не понимаю, что должна вычислять функция sg1. Значение сплайна в точке x, что ли? А зачем тогда у неё два аргумента? И зачем Вы пытаетесь сделать так, чтобы она возвращала массив?

Вы правы, n тут выступает в роли "масло-масленное".

Но по сути, убрав его, проблема не решается.

п.с. я делаю "подсетку" для того, чтобы сплайн можно было задать формулой

$f(x_j)\frac{x-x_{j+1}}{x_j-x_{j+1}}+f(x_{j+1} )\frac{x-x_j}{x_{j+1}-x_j}$

для j-того отрезка.

[_Ivana]

Сложности у вас не с методом а с его реализацией. Напишите корректно функцию сплайна. Потом вызовите её и исходную функцию для подсетки.

(Оффтоп)

Вот вам для сравнения и успокоения мой первый код на Матлабе - графики зависимости точности от шага сетки для разных сплайнов. А изначально я моделирование проводил в 1С (пока Матлаб не скачал) - хотите, дам листинг 1С кода?

function Main()
clear;

N_points_per = 600;
dt = 2*pi/N_points_per;

k_N = sqrt(sqrt(2));
N = 2;

num_of_algs = 7;
m_error = zeros(num_of_algs + 1, 1);

num_of_points = 42;
ip_0 = num_of_points/2;
y_i = zeros(1, num_of_points);
p_i = zeros(1, 22);
pp_i = zeros(1, 2);

count = 0;
while N < 129
   t_N = 2*pi/N;
   error_cur = zeros(num_of_algs + 1, 1);

   phase = 0;
   while phase < 2*pi
      t = 0;
      t_abs = 0;
      t_left_int = 0;
      for ii = 1:num_of_points y_i(ii) = sin(t_left_int + (ii - ip_0)*t_N + phase); end
      for ii = 1:22 p_i(ii) = Ivana_1_deriv(y_i, ii); end
      for ii = 1:2 pp_i(ii) = Ivana_1_deriv(p_i, ii); end

      while t_abs < 3*pi
         if t >= t_N
            t = t - t_N;
            t_left_int = t_left_int + t_N;
   
%            y_i(1) = [];
%            y_i(num_of_points) = sin(t_left_int + (num_of_points - ip_0)*t_N + phase);

            for ii = 1:(num_of_points - 1) y_i(ii) = y_i(ii + 1); end
            y_i(num_of_points) = sin(t_left_int + (num_of_points - ip_0)*t_N + phase);

            for ii = 1:21 p_i(ii) = p_i(ii + 1); end
            p_i(22) = Ivana_1_deriv(y_i, 22);

            pp_i(1) = pp_i(2);
            pp_i(2) = Ivana_1_deriv(p_i, 2);
         end

         y_sin  = sin(t_abs + phase);
         t_norm = t / t_N;

         y_int = y_i(0 + ip_0) + (y_i(1 + ip_0) - y_i(0 + ip_0))*t_norm;
         error_cur(2) = max(error_cur(2), abs(y_int - y_sin));

         a3 = (y_i(2 + ip_0) - y_i(-1 + ip_0))/6 + (y_i(0 + ip_0) - y_i(1 + ip_0))/2;
         a1 = (y_i(2 + ip_0) - y_i(0 + ip_0))/2 - a3;
         a2 = y_i(2 + ip_0) - y_i(1 + ip_0) - a1 - a3;
         a0 = y_i(1 + ip_0);
         t_norm_1 = t_norm - 1;
         y_int = a3*t_norm_1; y_int = (y_int + a2)*t_norm_1; y_int = (y_int + a1)*t_norm_1; y_int = y_int + a0;
         error_cur(3) = max(error_cur(3), abs(y_int - y_sin));

         p0 = (2*(y_i(1 + ip_0) - y_i(-1 + ip_0)) - (y_i(2 + ip_0) - y_i(-2 + ip_0))/4)/3;
         p1 = (2*(y_i(2 + ip_0) - y_i(0 + ip_0) ) - (y_i(3 + ip_0) - y_i(-1 + ip_0))/4)/3;
         a0 = y_i(0 + ip_0);
         a1 = p0;
         a3 = p0 + p1 + 2*(y_i(0 + ip_0) - y_i(1 + ip_0));
         a2 = y_i(1 + ip_0) - a3 - a1 - a0;
         y_int = a3*t_norm; y_int = (y_int + a2)*t_norm; y_int = (y_int + a1)*t_norm; y_int = y_int + a0;
         error_cur(4) = max(error_cur(4), abs(y_int - y_sin));

         p0 = (y_i(1 + ip_0) - y_i(-1 + ip_0))*3/4 - (y_i(2 + ip_0) - y_i(-2 + ip_0))*3/20 + (y_i(3 + ip_0) - y_i(-3 + ip_0))*1/60;
         p1 = (y_i(2 + ip_0) - y_i(0 + ip_0) )*3/4 - (y_i(3 + ip_0) - y_i(-1 + ip_0))*3/20 + (y_i(4 + ip_0) - y_i(-2 + ip_0))*1/60;
         a0 = y_i(0 + ip_0);
         a1 = p0;
         a3 = p0 + p1 + 2*(y_i(0 + ip_0) - y_i(1 + ip_0));
         a2 = y_i(1 + ip_0) - a3 - a1 - a0;
         y_int = a3*t_norm; y_int = (y_int + a2)*t_norm; y_int = (y_int + a1)*t_norm; y_int = y_int + a0;
         error_cur(5) = max(error_cur(5), abs(y_int - y_sin));

         p0 = 0   + (y_i(1 + ip_0) - y_i(-1 + ip_0))*5/6....
            - (y_i(2 + ip_0) - y_i(-2 + ip_0))*5/21....
            + (y_i(3 + ip_0) - y_i(-3 + ip_0))*5/84....
             - (y_i(4 + ip_0) - y_i(-4 + ip_0))*5/504....
            + (y_i(5 + ip_0) - y_i(-5 + ip_0))*1/1260;
         p1 = 0   + (y_i(2 + ip_0) - y_i( 0 + ip_0))*5/6....
            - (y_i(3 + ip_0) - y_i(-1 + ip_0))*5/21....
            + (y_i(4 + ip_0) - y_i(-2 + ip_0))*5/84....
             - (y_i(5 + ip_0) - y_i(-3 + ip_0))*5/504....
            + (y_i(6 + ip_0) - y_i(-4 + ip_0))*1/1260;
         pp0 = 2*(0....
            + ((y_i(1 + ip_0) + y_i(-1 + ip_0)) / 2 - y_i(0 + ip_0))*3/2....
            - ((y_i(2 + ip_0) + y_i(-2 + ip_0)) / 2 - y_i(0 + ip_0))*3/20....
            + ((y_i(3 + ip_0) + y_i(-3 + ip_0)) / 2 - y_i(0 + ip_0))*1/90);
         pp1 = 2*(0....
            + ((y_i(2 + ip_0) + y_i( 0 + ip_0)) / 2 - y_i(1 + ip_0))*3/2....
            - ((y_i(3 + ip_0) + y_i(-1 + ip_0)) / 2 - y_i(1 + ip_0))*3/20....
            + ((y_i(4 + ip_0) + y_i(-2 + ip_0)) / 2 - y_i(1 + ip_0))*1/90);
         y_int = Ivana_spline_E5_exact_deriv(t_norm, y_i(0 + ip_0), y_i(1 + ip_0), p0, p1, pp0, pp1);
         error_cur(6) = max(error_cur(6), abs(y_int - y_sin));

         y_int = Ivana_spline_E5_exact_deriv(t_norm, y_i(0 + ip_0), y_i(1 + ip_0), p_i(11), p_i(12), pp_i(1), pp_i(2));
         error_cur(7) = max(error_cur(7), abs(y_int - y_sin));

         p0 = t_N*cos(t_left_int  + phase);
         p1 = t_N*cos(t_left_int + t_N + phase);
         pp0 = -t_N*t_N*y_i(0 + ip_0);
         pp1 = -t_N*t_N*y_i(1 + ip_0);
         y_int = Ivana_spline_E5_exact_deriv(t_norm, y_i(0 + ip_0), y_i(1 + ip_0), p0, p1, pp0, pp1);
         error_cur(8) = max(error_cur(8), abs(y_int - y_sin));

         t = t + dt;
         t_abs = t_abs + dt;
      end

      phase = phase + 2*pi/13;
   end

   error_cur = log10(error_cur);
   error_cur(1) = log2(N);
   m_error = [m_error error_cur];   
      
   N = N*k_N;
   count = count + 1;
end

m_error(:, 1) = [];
plot(m_error(1, :), m_error(2, :), '-x', m_error(1, :), m_error(3, :), '-x', m_error(1, :), m_error(4, :), '-x', m_error(1, :), m_error(5, :), '-x'....
, m_error(1, :), m_error(6, :), '-x', m_error(1, :), m_error(7, :), '-x', m_error(1, :), m_error(8, :), '-x')
grid on
xlabel('логарифм по основанию 2 количества точек на период')
ylabel('логарифм максимальной абсолютной ошибки интерполяции')
title('Графики ошибки интерполяции от относительной частоты дискретизации')
legend('Лагранж-1', 'Лагранж-3', '[1/5, 1/5]', '[1/7, 1/7]', '[1/11/7, 1/11/7]', '[1/21/21, 1/21/21]', 'Эрмит-точные 1 и 2')

end

function r = Ivana_spline_E5_exact_deriv(t, y0, y1, p0, p1, pp0, pp1)
%
         a0 = y0;
         a1 = p0;
         a2 = pp0/2;
         k1 = y1 - a0 - a1 - a2;
         k2 = p1 - p0 - pp0;
         k3 = (pp1 - pp0)/2;
         a5 = k3 - 3*k2 + 6*k1;
         a4 = 7*k2 - 2*k3 - 15*k1;
         a3 = k1 - a5 - a4;
         r = a5*t; r = (r + a4)*t; r = (r + a3)*t; r = (r + a2)*t; r = (r + a1)*t; r = r + a0;
end

function r = Ivana_1_deriv(a, n)
%
i0 = n + 10;

r =  (a(i0 + 1) - a(i0 - 1))*10/11....
    -(a(i0 + 2) - a(i0 - 2))*15/44....
    +(a(i0 + 3) - a(i0 - 3))*20/143....
    -(a(i0 + 4) - a(i0 - 4))*15/286....
    +(a(i0 + 5) - a(i0 - 5))*12/715....
    -(a(i0 + 6) - a(i0 - 6))*15/3432....
    +(a(i0 + 7) - a(i0 - 7))*15/17017....
    -(a(i0 + 8) - a(i0 - 8))*5/38896....
    +(a(i0 + 9) - a(i0 - 9))*5/415701....
    -(a(i0 + 10) - a(i0 - 10))*1/1847560;

end

Сейчас я бы конечно написал пооптимальнее. Но все равно желающие могут скопипастить и посмотреть результат

[Alexeybk5]

_Ivana в сообщении #638953 писал(а):

Сложности у вас не с методом а с его реализацией. Напишите корректно функцию сплайна. Потом вызовите её и исходную функцию для подсетки.

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

a = -10.;
b = 10.;
n = 10.;
f[x_] := Exp[x]

h = (b - a)/n;
X[n_] := Table[{a + h *k, f[a + h*k]}, {k, 0, n}]
X[n] // MatrixForm

Функция сплайна и его значение в точке x=6

sg1[x_] := Module[{n, j, X = X[n]},
  For[j = 1, j <= n + 1, j++,
   If[x < (j)*h && x > (j - 1)*h, 
    s = Table[
      f[X[[j, 2]]]*(x - X[[j + 1, 1]])/(X[[j, 1]] - X[[j + 1, 1]]) + 
       f[X[[j + 1, 
           2]]]*((x - X[[j, 1]])/(X[[j + 1]] - X[[j, 1]])), {j, 1, 
       n}]]]
  ]

sg1[6]
s
то почему на выходе я получаю такой странный вектор, вместо просто числа ?
{5, {1.71388, -1.53769}}

p.ps
Да... Программа у Вас была короткая ))

[ИСН]

Пока Вы не уберёте этот бредовый цикл, я не скажу, в какой строчке программы у Вас из-за пропущенного индекса происходит вычитание между вектором и числом, которое Mathematica, тварь такая, слава её создателям, молчаливо интерпретирует как поэлементное вычитание, дающее тоже вектор - и дальше по всем операциям гордо несёт этот вектор.

[Alexeybk5]

ИСН в сообщении #638961 писал(а):

Пока Вы не уберёте этот бредовый цикл, я не скажу, в какой строчке программы у Вас из-за пропущенного индекса происходит вычитание между вектором и числом, которое Mathematica, тварь такая, слава её создателям, молчаливо интерпретирует как поэлементное вычитание, дающее тоже вектор - и дальше по всем операциям гордо несёт этот вектор.

А чем он плох ? Вроде бы обычный цикл.

Правда теперь вообще в s ничего не записывается.... Странная эта mathematica...

a = -10.;
b = 10.;
n = 10.;
f[x_] := Exp[x]

h = (b - a)/n;
X[n_] := Table[{a + h *k, f[a + h*k]}, {k, 0, n}]
X[n] // MatrixForm

sg1[x0_] := Module[{, y = x0, n, j, X = X[n]},
  For[j = 1, j <= n + 1, j++,
   If[y < (j)*h && y > (j - 1)*h, 
    s = f[X[[j, 2]]]*(y - X[[j + 1, 1]])/(X[[j, 1]] - X[[j + 1, 1]]) +
       f[X[[j + 1, 2]]]*((y - X[[j, 1]])/(X[[j + 1]] - X[[j, 1]]))]]
  ]

sg1[40]
s

[_Ivana]

Alexeybk5 в сообщении #638954 писал(а):

p.ps
Да... Программа у Вас была короткая ))

Зато рабочая и без глюков И делает что надо - строит сравнительные графики точности аппроксимации синуса от шага сетки семью различными сплайнами - от вашего же линейного до других поинтереснее, которые считаются не так тривиально.

ЗЫ не можете написать в матричном виде - напишите в скалярном через циклы, как и сделал я в своем коде.

ЗЗЫ тут была фраза про то, что мне не понравилось одно место в вашем коде, но вдумчиво прочитав сообщение ИСН, я её стер. Не скажем, пока не уберете бредовый цикл

[ИСН]

Один студент © тоже выучил циклы и думал, что больше ничего не надо. Так и писал. Скажем, чтобы поделить с остатком 1000 на 17, он устраивал цикл:

For[i = 1, i <= 1000, i++,
  If[(i)*17 < 1000 && (i + 1)*17 > 1000,
   s = i
   ]
  ];
s

Потом он попал в армию. Там ему дали лом и велели копать яму.
- Но так же неудобно! - возмутился бывший студент.
- Копай, ..... ............., а то ..... ..... - сказали ему.

[Alexeybk5]

Согласен, цикл в мусорку.

Я постарался изменить код , но теперь вообще не работает ...

a = -10.;
b = 10.;
n = 10.;
f[x_] := Exp[x]
h = (b - a)/n;
X[n_] := Table[{a + h*k, f[a + h*k]}, {k, 0, n}];

sgg[x_] := 
 Block[{X = X[n]}, 
  f[X[[j, 2]]] (y - X[[j + 1, 1]])/(X[[j, 1]] - X[[j + 1, 1]]) + 
   f[X[[j + 1, 2]]] ((y - X[[j, 1]])/(X[[j + 1, 1]] - X[[j, 1]]))]

apprg[x_] := 
 Block[{$MaxExtraPrecision = Infinity, y = SetPrecision[x, 500]}, 
  N[sgg[y] // Expand, 500]]

Plot[{apprg[x], f[x]}, {x, a, b}, ImageSize -> 600, 
 PlotStyle -> {Dashed, Dotted, Red, Thick}]

data[n_] := Table[{x, Abs[apprg[x] - f[x]]}, {x, a, b, 1/100}]

ListLinePlot[Table[Tooltip[data[n], n], {n, 1}], PlotRange -> All, 
 ImageSize -> 600, PlotStyle -> Thick]

[ИСН]

Уберите все сложные конструкции. Уберите нахрен Block, вычисляйте массив X один раз в начале. И больше ничего сложного не трогайте, пока не перестанете попадать вилкой в глаз. Вот это f[X[[j, 2]]] - это что такое, например?

[Alexeybk5]

Получается, что это функция от функции... Действительно, X[[j, 2]] уже и так значене функции в узле))

Ещё одна попытка

a = -10.;
b = 10.;
n = 10.;
f[x_] := Exp[x]
h = (b - a)/n;
X[n_] := Table[{a + h*k, f[a + h*k]}, {k, 0, n}];

X = X[n]
m=10.;
ss := Table[{a + (h/m)*k, f[a + (h/m)*k]}, {k, 0, n*m}]
ss
sp1 := Table[ss[[i, 1]], 
  X[[j, 2]] (ss[[i, 1]] - X[[j + 1, 1]])/(X[[j, 1]] - X[[j + 1, 1]]) + 
   X[[j + 1, 
      2]] ((ss[[j, 1]] - X[[j, 1]])/(X[[j + 1, 1]] - X[[j, 1]])), {i, 
   0, n, h/m}, {j, 0, n}]
sp1

Ну или альтернатива Table
ещё один кривой, двойной, цикл

For[j = 1, j <= n, j++
  For[i = j, i <= m*j, i++ ,
   
   spsp = 
    X[[j, 2]] (ss[[i, 1]] - X[[j + 1, 1]])/(X[[j, 1]] - 
         X[[j + 1, 1]]) + 
     X[[j + 1, 
        2]] ((ss[[j, 1]] - X[[j, 1]])/(X[[j + 1, 1]] - X[[j, 1]]))]]
spsp

Но он то почему не работает ?
Ведь первый For присвает j=1
после чего второй For должен проводить операцию для i= 1 до i=1*m=1*10=10
после чего заносить 10 строк в вектор spsp
и опять запускать первый For для j = 2 и т.д. Но ведь сейчас же верно записан, почему не работает...?

[ИСН]

Уберите второй Table (почему в нём двойной цикл?), и форов тоже не треба. Оформите сплайн в виде одной функции, которая берёт один долбаный аргумент и возвращает одно #$% значение. Один и одно.

[Alexeybk5]

Ну а как же потом тогда разницу апроксимации функции сплайном определить ?
Ведь надо , положим, функцию экспоненты апроксимировать не одним сплайном 1 порядка, а n сплайнами 1 порядка.
....
Т.е. Вы предлагаете, строить одну сетку, на которой находить значения функции, после чего с помощью формулы найти значения в этих узлах сплайна
$f(x_j)\frac{x-x_{j+1}}{x_j-x_{j+1}}+f(x_{j+1} )\frac{x-x_j}{x_{j+1}-x_j}$
И какой же $x$ сюда подставлять в таком случае? Ведь мы просто получим в итоге ровно такую же таблицу, что и таблица значений функции в этих точках... Т.е. ничего не получим...
Всё равно надо строить подсетку, а для этого и нужен пторой итератор i...
Но опять же, а как потом найти разницу апроксимации ?

[ИСН]

Да, надо строить подсетку. Но Вы не можете её строить, пока не имеете функции, которая будет по произвольному числу x возвращать значение сплайна в этой точке. Сделайте эту функцию. Сделайте её. Оценок не надо, оценки потом.

[Alexeybk5]

ИСН в сообщении #639152 писал(а):

Да, надо строить подсетку. Но Вы не можете её строить, пока не имеете функции, которая будет по произвольному числу x возвращать значение сплайна в этой точке. Сделайте эту функцию. Сделайте её. Оценок не надо, оценки потом.

а если так сделать ?

a = -10.;
b = 10.;
n = 10.;
f[x_] := Exp[x]
h = (b - a)/n;
X[n_] := Table[{a + h*k, f[a + h*k]}, {k, 0, n}];
X = X[n]
m = 10.;
ss := Table[{a + (h/m)*k, f[a + (h/m)*k]}, {k, 0, n*m}]
ss
sbl = X[[i, 
     2]] (ss[[i, 1]] - X[[i*m + 1, 1]])/(X[[i*m, 1]] - 
      X[[i*m + 1, 1]]) + 
  X[[i*m + 1, 
     2]] ((ss[[i*m, 1]] - X[[i*m, 1]])/(X[[i*m + 1, 1]] - 
       X[[i*m, 1]])), {i, 0, n, h/m}

Просто функция sbl в которую должно вноситься реузльтат работы с двумя таблицами.... Ну а как иначе-то реализовать....

п.с. функция у меня есть, но вы её забраковали

a = -10.;
b = 10.;
n = 10.;
f[x_] := Exp[x]
h = (b - a)/n;
X[n_] := Table[{a + h*k, f[a + h*k]}, {k, 0, n}];

sg1[x0_] :=
Module[{y = x0, j, X = X[n]},
For[j = 1, j <= n + 1, j++,
If[(y < j h) && (y > (j - 1) h),
s = 
f[X[[j, 2]]] (y - X[[j + 1, 1]])/(X[[j, 1]] - X[[j + 1, 1]]) +
f[X[[j + 1, 
2]]] ((y - X[[j, 1]])/(X[[j + 1, 1]] - X[[j, 1]]))
];
];
];
sg1[4.4]
s

m = 10.;
sss := Table[{a + (h/m)*k}, {k, 0, n*m}]
sss
sg1[sss]
s

Положим, что я готов ломом копать яму., как сделать это до конца? Ведь научившись копать ломом, потом лопатой я ж точно смогу.

[ИСН]

Вашу функцию я забраковал потому, что во-первых, она возвращает бред, а во-вторых, у неё внутри этот чёртов цикл, где Вы находите результат деления - подбором. Если не это копание ломом, то что?