Вычисление функций рядом Маклорена

% Вычисление функций рядом Маклорена
% d-derivatives (производные функции вычисленные в нуле)
% x- аргумент вычисляемой функции
% Например:
%     d=ones(32,1);
%     x=1;
%     y=fmacloren01(d,x);
% 
function y=fmacloren01(d,x)
y=0;
p=1;
f=1;
n=length(d);
if n<1 
    return;
end;
for k=1:n
    y=y+(d(k)/f)*p;
    f=f*k;
    p=p*x;
end;


% Вычисление функции SIN(x) рядом Маклорена
% x- аргумент вычисляемой функции
% Например:
%     x=pi/2;
%     y=fsin01(x);
%
function y=fsin01(x)
y=0;  % Предустановка результата
p=1;  % x - в очередной степени
f=1;  % Очередное значение факториала
% d-derivatives (производные функции вычисленные в нуле)
d=[0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1];
n=length(d);
for k=1:n
    y=y+(d(k)/f)*p;
    f=f*k;
    p=p*x;
end;


% Вычисление функции SIN(x) рядом Маклорена
% (Оптимизация алгоритма для исключения вычислений с нулевыми производными)
% x- аргумент вычисляемой функции
% Например:
%     x=pi/2;
%     y=fsin02(x);
% 
function y=fsin02(x)
y=0;  % Предустановка результата
p=1;  % x - в очередной степени
f=1;  % Очередное значение факториала
% d-derivatives (производные функции вычисленные в нуле)
d=[0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1];
n=length(d);
for k=1:n
    if d(k) ~= 0
       y=y+(d(k)/f)*p;
    end;   
    f=f*k;
    p=p*x;    
end;


% Вычисление функции SIN(x) рядом Маклорена
% (Оптимизация алгоритма с ограничением длины ряда по достижению ошибки delta)
% x- аргумент вычисляемой функции
% Например:
%     x=pi/2;
%     y=fsin03(x,1e-10);
% 
function y=fsin03(x,delta)
y=0;  % Предустановка результата
p=1;  % x - в очередной степени
f=1;  % Очередное значение факториала
% d-derivatives (производные функции вычисленные в нуле)
d1=[0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1];
d2=[0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1];
d=[d1,d2];
s=realmax;
n=length(d);
for k=1:n
    if d(k)~=0
       y=y+(d(k)/f)*p;
       if abs(s-y)<delta
          disp(strcat('Длина ряда k=', int2str(k), ', макс. длина n=', int2str(n)));    
          return;
       end;
       s=y;
    end;   
    f=f*k;
    p=p*x;    
end;


% Вычисление функции COS(x) рядом Маклорена
% (Оптимизация алгоритма для исключения вычислений с нулевыми производными)
% x- аргумент вычисляемой функции
% Например:
%     x=pi/2;
%     y=fcos02(x);
% 
function y=fcos02(x)
y=0;  % Предустановка результата
p=1;  % x - в очередной степени
f=1;  % Очередное значение факториала
% d-derivatives (производные функции вычисленные в нуле)
d1=[1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 ];
d2=[1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 ];
d=[d1,d2];
n=length(d);
for k=1:n
    if d(k)~=0
       y=y+(d(k)/f)*p;  
    end;   
    f=f*k;
    p=p*x;    
end;



% Вычисление функции COS(x) рядом Маклорена
% (Оптимизация алгоритма с ограничением длины ряда по достижению ошибки delta)
% x- аргумент вычисляемой функции
% Например:
%     x=pi/2;
%     y=fcos03(x,1e-10);
%
function y=fcos03(x,delta)
y=0;  % Предустановка результата
p=1;  % x - в очередной степени
f=1;  % Очередное значение факториала
% d-derivatives (производные функции вычисленные в нуле)
d1=[1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 ];
d2=[1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 ];
d=[d1,d2];
s=inf;
n=length(d);
for k=1:n
    if d(k)~=0
       y=y+(d(k)/f)*p;
       if abs(s-y)<delta
          disp(strcat('Длина ряда k=', int2str(k), ', макс. длина n=', int2str(n)));   
          return;
       end;
       s=y;
    end;   
    f=f*k;
    p=p*x;    
end;


% Вычисление функции COS(x) рядом Маклорена
% (Расширенный диапазон аргумента)
% x- аргумент вычисляемой функции
% Например:
%     x=pi/2;
%     y=fcos04(x);
% 
function y=fcos04(x)
y=0;  % Предустановка результата
p=1;  % x - в очередной степени
f=1;  % Очередное значение факториала
% Вычислим отаток от деления на период
x1=rem(x,2*pi);
% d-derivatives (производные функции вычисленные в нуле)
d1=[1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 ];
d2=[1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 ];
d=[d1,d2];
n=length(d);
for k=1:n
    if d(k)~=0
       y=y+(d(k)/f)*p;  
    end;   
    f=f*k;
    p=p*x1;    
end;



% Тестирование точности функции fcos02(x) по библиотечной cos(x)
% Синтаксис: [delta]=ftestcos02(xb, xe)
% xb - начало диапазона аргумента
% xe - конец диапазона аргумента
% Пример вызова: delta=ftestcos02(0, pi);
% 
function [delta]=ftestcos02(xb, xe)
x=[xb:pi/100:xe];
n=length(x);
for k=1:n
   y1(k)= fcos02(x(k));  
   y2(k)= cos(x(k));
   delta(k)= y1(k)-y2(k);
end;
grid on;
plot(x,delta);


% Тестирование точности функции fcos04(x) по библиотечной cos(x)
% (Расширенный диапазон аргумента)
% Синтаксис: [delta]=ftestcos04(xb, xe)
% xb - начало диапазона аргумента
% xe - конец диапазона аргумента
% Пример вызова: delta=ftestcos04(-8*pi, 8*pi);
% 
function [delta]=ftestcos04(xb, xe)
x=[xb:pi/100:xe];
n=length(x);
for k=1:n
   y1(k)= fcos04(x(k));  
   y2(k)= cos(x(k));
   delta(k)= y1(k)-y2(k);
end;
grid on;
plot(x,delta);


% Вычисление производных табличной SIN
% Вычисляются следующие производные:
% Первая    - (blue график)  =   cos(x)
% Вторая    - (green график) = - sin(x)
% Третья    - (Red график)   = - cos(x)
% Четвертая - (Cyan график)  =   sin(x)
%
function deriv_sin01
x=[0:pi/500:2*pi];
y=sin(x);
grid on;
hold on;

for k=1:(length(x)-1)
    z1(k)=(y(k+1)-y(k))/(x(k+1)-x(k));
    a(k)=x(k);
end; 
plot(a,z1);
clear a;
for k=1:(length(x)-2)
    z2(k)=(z1(k+1)-z1(k))/(x(k+1)-x(k));
    a(k)=x(k);
end;
plot(a,z2,'g');
clear a;
for k=1:(length(x)-3)
    z3(k)=(z2(k+1)-z2(k))/(x(k+1)-x(k));
    a(k)=x(k);
end;
plot(a,z3,'r');
clear a;
for k=1:(length(x)-4)
    z4(k)=(z3(k+1)-z3(k))/(x(k+1)-x(k));
    a(k)=x(k);
end;
plot(a,z4,'c');


% Вычисление производных зашумленной табличной SIN 
% Синтаксис: deriv_sin2_noise(noisescale)
% Например:
%            deriv_sin01_noise(1.e-4);
% 
% Вычисляются следующие производные:
% Первая    - (blue график)  =   cos(x)
% Вторая    - (green график) = - sin(x)
%
function deriv_sin2_noise(noisescale)
x=[0:pi/500:2*pi];
y=sin(x);
% Сформировать и отмасштабировать шум
noise=rand(length(y),1);
noise=noise.*noisescale;
noise=noise';
disp('Макс. амплитуда сигнала = 1');
ampnoise = max(abs(noise));
disp(strcat('Макс. амплитуда шума = ', num2str(ampnoise)));
% Сформировать зашумленную табличную функцию
y=y+noise;
clf;
grid on;
hold on;
for k=1:(length(x)-1)
    z1(k)=(y(k+1)-y(k))/(x(k+1)-x(k));
    a(k)=x(k);
end; 
plot(a,z1);
clear a;
for k=1:(length(x)-2)
    z2(k)=(z1(k+1)-z1(k))/(x(k+1)-x(k));
    a(k)=x(k);
end;
plot(a,z2,'g');


% Вычисление производных зашумленной табличной SIN 
% Синтаксис: deriv_sin01_noise(noisescale)
% Например:
%            deriv_sin3_noise(1.e-7);
% 
% Вычисляются следующие производные:
% Первая    - (blue график)  =   cos(x)
% Вторая    - (green график) = - sin(x)
% Третья    - (Red график)   = - cos(x)
%
function deriv_sin3_noise(noisescale)
x=[0:pi/500:2*pi];
y=sin(x);
% Сформировать и отмасштабировать шум
noise=rand(length(y),1);
noise=noise.*noisescale;
noise=noise';
disp('Макс. амплитуда сигнала = 1');
ampnoise = max(abs(noise));
disp(strcat('Макс. амплитуда шума = ', num2str(ampnoise)));
% Сформировать зашумленную табличную функцию
y=y+noise;
clf;
grid on;
hold on;
for k=1:(length(x)-1)
    z1(k)=(y(k+1)-y(k))/(x(k+1)-x(k));
    a(k)=x(k);
end; 
plot(a,z1);
clear a;
for k=1:(length(x)-2)
    z2(k)=(z1(k+1)-z1(k))/(x(k+1)-x(k));
    a(k)=x(k);
end;
plot(a,z2,'g');
clear a;
for k=1:(length(x)-3)
    z3(k)=(z2(k+1)-z2(k))/(x(k+1)-x(k));
    a(k)=x(k);
end;
plot(a,z3,'r');


% Вычисление производных зашумленной табличной SIN 
% Синтаксис: deriv_sin01_noise(noisescale)
% Например:
%            deriv_sin4_noise(1.e-10);
% 
% Вычисляются следующие производные:
% Первая    - (blue график)  =   cos(x)
% Вторая    - (green график) = - sin(x)
% Третья    - (Red график)   = - cos(x)
% Четвертая - (Cyan график)  =   sin(x)
%
function deriv_sin4_noise(noisescale)
x=[0:pi/500:2*pi];
y=sin(x);
% Сформировать и отмасштабировать шум
noise=rand(length(y),1);
noise=noise.*noisescale;
noise=noise';
disp('Макс. амплитуда сигнала = 1');
ampnoise = max(abs(noise));
disp(strcat('Макс. амплитуда шума = ', num2str(ampnoise)));
% Сформировать зашумленную табличную функцию
y=y+noise;
clf;
grid on;
hold on;
for k=1:(length(x)-1)
    z1(k)=(y(k+1)-y(k))/(x(k+1)-x(k));
    a(k)=x(k);
end; 
plot(a,z1);
clear a;
for k=1:(length(x)-2)
    z2(k)=(z1(k+1)-z1(k))/(x(k+1)-x(k));
    a(k)=x(k);
end;
plot(a,z2,'g');
clear a;
for k=1:(length(x)-3)
    z3(k)=(z2(k+1)-z2(k))/(x(k+1)-x(k));
    a(k)=x(k);
end;
plot(a,z3,'r');
clear a;
for k=1:(length(x)-4)
    z4(k)=(z3(k+1)-z3(k))/(x(k+1)-x(k));
    a(k)=x(k);
end;
plot(a,z4,'c');