Первый код на Матлабе :), желающие ругают, я работаю над ошибками Опции

[_Ivana]

Собственно, делюсь радостью - я наконец-то обрел и поставил себе Матлаб 2012а Под это дело купил себе на днях новый комп даже Скачал несколько обучающих пдф-ок, прочитал пока только одну - Get start и сразу загорелся идеей смоделировать мои предыдущие эксперименты в 1С по расчету максимальной ошибки различных методов интерполяции - даже получилось. Мыслить в категориях массивов и матриц пока не привык, поэтому в коде скорее всего много глупостей и неоптимальностей, лишних переменных и т.п., не получается вывести график без линеаризации масштаба по оси Х (он сам его линеаризует) и т.д. Буду играться дальше, а пока желающим представляю сам код - для замечаний и предложений, за которые заранее спасибо

CODE

clear;

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

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

m_N = 0;
m_count = 0;
count = 0;
m_error_L1 = 0;
m_error_L3 = 0;
m_error_E5 = 0;
m_error_E7 = 0;

while N < 259
t_N = 2*pi/N;

error_L1 = 0;
error_L3 = 0;
error_E5 = 0;
error_E7 = 0;

phase = 0;
while phase < 2*pi

t = 0;
t_abs = 0;

t_left_int = 0;
y_3 = sin(t_left_int - 3*t_N + phase);
y_2 = sin(t_left_int - 2*t_N + phase);
y_1 = sin(t_left_int - t_N + phase);
y0 = sin(t_left_int + phase);
y1 = sin(t_left_int + t_N + phase);
y2 = sin(t_left_int + 2*t_N + phase);
y3 = sin(t_left_int + 3*t_N + phase);
y4 = sin(t_left_int + 4*t_N + phase);

% x = 0;
% y_sin = 0;
% y_int = 0;

while t_abs < 3*pi
if t >= t_N
t = t - t_N;
t_left_int = t_left_int + t_N;

y_3 = y_2;
y_2 = y_1;
y_1 = y0;
y0 = y1;
y1 = y2;
y2 = y3;
y3 = y4;
y4 = sin(t_left_int + 4*t_N + phase);
end

% x = [x t_abs];
% y_sin = [y_sin sin(t_abs + phase)];
y_sin = sin(t_abs + phase);
t_norm = t / t_N;

y_int_cur = y0 + (y1 - y0)*t_norm;
error_L1_cur = abs(y_int_cur - y_sin); if error_L1_cur > error_L1 error_L1 = error_L1_cur; end

a3 = (y2 - y_1)/6 + (y0 - y1)/2;
a1 = (y2 - y0)/2 - a3;
a2 = y2 - y1 - a1 - a3;
a0 = y1;
t_norm_1 = t_norm - 1;
y_int_cur = a3*t_norm_1; y_int_cur = (y_int_cur + a2)*t_norm_1; y_int_cur = (y_int_cur + a1)*t_norm_1; y_int_cur = y_int_cur + a0;
% y_int = [y_int y_int_cur];
error_L3_cur = abs(y_int_cur - y_sin); if error_L3_cur > error_L3 error_L3 = error_L3_cur; end

p0 = (2*(y1-y_1) - (y2 - y_2)/4)/3;
p1 = (2*(y2-y0) - (y3 - y_1)/4)/3;
a0 = y0;
a1 = p0;
a3 = p0 + p1 + 2*(y0 - y1);
a2 = y1 - a3 - a1 - a0;
y_int_cur = a3*t_norm; y_int_cur = (y_int_cur + a2)*t_norm; y_int_cur = (y_int_cur + a1)*t_norm; y_int_cur = y_int_cur + a0;
error_E5_cur = abs(y_int_cur - y_sin); if error_E5_cur > error_E5 error_E5 = error_E5_cur; end

p0 = (y1-y_1)*3/4 - (y2 - y_2)*3/20 + (y3 - y_3)/60;
p1 = (y2-y0 )*3/4 - (y3 - y_1)*3/20 + (y4 - y_2)/60;
a0 = y0;
a1 = p0;
a3 = p0 + p1 + 2*(y0 - y1);
a2 = y1 - a3 - a1 - a0;
y_int_cur = a3*t_norm; y_int_cur = (y_int_cur + a2)*t_norm; y_int_cur = (y_int_cur + a1)*t_norm; y_int_cur = y_int_cur + a0;
error_E7_cur = abs(y_int_cur - y_sin); if error_E7_cur > error_E7 error_E7 = error_E7_cur; end

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

% x(1) = [];
% y_sin(1) = [];
% y_int(1) = [];
% plot(x, y_sin, x, y_int)

phase = phase + 2*pi/13;
end

m_N = [m_N N];
m_count = [m_count count];
m_error_L1 = [m_error_L1 log10(error_L1)];
m_error_L3 = [m_error_L3 log10(error_L3)];
m_error_E5 = [m_error_E5 log10(error_E5)];
m_error_E7 = [m_error_E7 log10(error_E7)];

N = N*k_N;
count = count + 1;
end

m_N(1) = [];
m_count(1) = [];
m_error_L1(1) = [];
m_error_L3(1) = [];
m_error_E5(1) = [];
m_error_E7(1) = [];
plot(m_N, m_error_L1, '-o', m_N, m_error_L3, '-o', m_N, m_error_E5, '-o', m_N, m_error_E7, '-o')
grid on
figure
plot(m_count, m_error_L1, '-o', m_count, m_error_L3, '-o', m_count, m_error_E5, '-o', m_count, m_error_E7, '-o')
grid on

[ViKo]

для замечаний и предложений, за которые заранее спасибо

По графикам не понять, где что. Их же, вроде подписывать можно.

[_Ivana]

Ну да, раз уж взялся выводить графики, то надо делать это красиво - с заголовками, подписями осей и легендой... Но с этим то как раз сложностей не должно возникнуть - команды для всего вышеперечисленного существуют и известны. С этим графиком меня больше волнует как сделать логарифмический масштаб по оси Х или же просто подписать её из отдельного массива не равноотстоящих отсчетов при равноотстоящем выводе - просто по порядковому номеру значения (для этого 2 графика для сравнения и вывожу).
ЗЫ сейчас уже доработал модель, вписал ещё несколько сплайнов, сетку исходных точек, итоговые результаты и некоторые промежуточные данные считаю в векторах и матрицах.
ЗЗЫ смайлики пришлось отключить, а то улыбчивый движок форума их и в теге CodeBox выводил - а там у меня двоеточий со скобками хватает :)

CODE

clear;

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

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

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

num_of_points = 12;
ip_0 = num_of_points/2;
y_i = zeros(1, num_of_points);

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

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);
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);
a0 = y_i(0 + ip_0);
a1 = p0;
a2 = pp0/2;
k1 = y_i(1 + ip_0) - 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;
y_int = a5*t_norm; y_int = (y_int + a4)*t_norm; y_int = (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(6) = max(error_cur(6), 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);
a0 = y_i(0 + ip_0);
a1 = p0;
a2 = pp0/2;
k1 = y_i(1 + ip_0) - 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;
y_int = a5*t_norm; y_int = (y_int + a4)*t_norm; y_int = (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(7) = max(error_cur(7), 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) = count;
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')
grid on
xlabel('количество точек на период, N (от 2 до 128, числам не верить :))')
ylabel('логарифм максимальной абсолютной ошибки интерполяции')
title('Графики ошибки интерполяции от относительной частоты дискретизации')
legend('Лагранж-1', 'Лагранж-3', '[1/5, 1/5]', '[1/7, 1/7]', '[1/11/7, 1/11/7]', 'Эрмит-точные 1 и 2')