«Оптимал квадратур формулалар орқали Фурье интегралларини тақрибий ҳисоблаш алгоритми» дастурий тизими

'),w2,real(F0_2'))
legend('FT', 'OQF (h=0.1)', 'OQF (h=0.01)', 'OQF (h=0.001)')
axis([a b - 5 5])
title('$F[f_0]$', 'Interpreter', 'latex');

subplot(3, 2, 4)

plot(w1, abs(real(F0_1
')-real(E0_1)))
title('$|\Re R_{f_{0,\omega}}[-1,1]| (h=0.1)$ ', 'Interpreter', 'latex');

subplot(3, 2, 5)

plot(w, abs(real(F0
')-real(E0)))
title('$|\Re R_{f_{0,\omega}}[-1,1]| (h=0.01)$ ', 'Interpreter', 'latex');

subplot(3, 2, 6)

plot(w2, abs(real(F0_2
')-real(E0_2)))
title('$|\Re R_{f_{0,\omega}}[-1,1]| (h=0.001)$ ', 'Interpreter', 'latex');

## OQF and FFT[fsin]
figure
subplot(3, 2, [1, 2])
plot(xx, fsin)
axis([-1.5 1.5 - 1.5 1.5])
title('$f_{\sin}$', 'Interpreter', 'latex');

subplot(3, 2, 3)

plot(w, real(E1), w1, real(F1_1), w, real(F1
'),w2,real(F1_2'))
legend('FT', 'OQF (h=0.1)', 'OQF (h=0.01)', 'OQF (h=0.001)')
axis([a b - 5 5])
title('$F[f_{\sin}]$', 'Interpreter', 'latex');

subplot(3, 2, 4)

plot(w1, abs(real(F1_1
')-real(E1_1)))
title('$|\Re R_{f_{{\sin},\omega}}[-1,1]| (h=0.1)$ ', 'Interpreter', 'latex');

subplot(3, 2, 5)

plot(w, abs(real(F1
')-real(E1)))
title('$|\Re R_{f_{{\sin},\omega}}[-1,1]| (h=0.01)$ ', 'Interpreter', 'latex');

subplot(3, 2, 6)

plot(w2, abs(real(F1_2
')-real(E1_2)))
title('$|\Re R_{f_{{\sin},\omega}}[-1,1]| (h=0.001)$ ', 'Interpreter', 'latex');

# OQF and FFT[fexp]
figure
subplot(3, 2, [1, 2])
plot(xx, fexp)
axis([-1.5 1.5 - 1.5 1.5])
title('$f_{\exp}$', 'Interpreter', 'latex');

subplot(3, 2, 3)

plot(w, real(E2), w1, real(F2_1), w, real(F2
'),w2,real(F2_2'))
legend('FT', 'OQF (h=0.1)', 'OQF (h=0.01)', 'OQF (h=0.001)')
axis([a b - 5 5])
title('$F[f_{\exp}]$', 'Interpreter', 'latex');

subplot(3, 2, 4)

plot(w1, abs(real(F2_1
')-real(E2_1)))
title('$|\Re R_{f_{{\exp},\omega}}[-1,1]| (h=0.1)$ ', 'Interpreter', 'latex');

subplot(3, 2, 5)

plot(w, abs(real(F2
')-real(E2)))
title('$|\Re R_{f_{{\exp},\omega}}[-1,1]| (h=0.01)$ ', 'Interpreter', 'latex');

subplot(3, 2, 6)

plot(w2, abs(real(F2_2
')-real(E2_2)))
title('$|\Re R_{f_{{\exp},\omega}}[-1,1]| (h=0.001)$ ', 'Interpreter', 'latex');

# OQF and FFT[flinear]
figure
subplot(3, 2, [1, 2])
plot(xx, flinear)
axis([-1.5 1.5 - 1.5 1.5])
title('$f_{linear}$', 'Interpreter', 'latex');

subplot(3, 2, 3)

plot(w, real(E11), w1, real(F11_1), w, real(F11
'),w2,real(F11_2'))
legend('FT', 'OQF (h=0.1)', 'OQF (h=0.01)', 'OQF (h=0.001)')
axis([a b - 5 5])
title('$F[f_linear]$', 'Interpreter', 'latex');

subplot(3, 2, 4)

plot(w1, abs(real(F11_1
')-real(E11_1)))
title('$|\Re R_{f_linear},\omega}}[-1,1]| (h=0.1)$ ', 'Interpreter', 'latex');

subplot(3, 2, 5)

plot(w, abs(real(F11
')-real(E11)))
title('$|\Re R_{f_{linear,\omega}}[-1,1]| (h=0.01)$ ', 'Interpreter', 'latex');

subplot(3, 2, 6)

plot(w2, abs(real(F11_2
')-real(E11_2)))
title('$|\Re R_{f_linear},\omega}}[-1,1]| (h=0.001)$ ', 'Interpreter', 'latex');

# OQF and FFT[fquad]
figure
subplot(3, 2, [1, 2])
plot(xx, fquad)
axis([-1.5 1.5 - 1.5 1.5])
title('$f_{\quad}$', 'Interpreter', 'latex');

subplot(3, 2, 3)

plot(w, real(E21), w1, real(F21_1), w, real(F21
'),w2,real(F21_2'))
legend('FT', 'OQF (h=0.1)', 'OQF (h=0.01)', 'OQF (h=0.001)')
axis([a b - 5 5])
title('$F[f_{\quad}]$', 'Interpreter', 'latex');

subplot(3, 2, 4)

plot(w1, abs(real(F21_1
')-real(E21_1)))
title('$|\Re R_{f_{{\quad},\omega}}[-1,1]| (h=0.1)$ ', 'Interpreter', 'latex');

subplot(3, 2, 5)

plot(w, abs(real(F21
')-real(E21)))
title('$|\Re R_{f_{{\exp},\omega}}[-1,1]| (h=0.01)$ ', 'Interpreter', 'latex');

subplot(3, 2, 6)

plot(w2, abs(real(F21_2
')-real(E21_2)))
title('$|\Re R_{f_{{\quad},\omega}}[-1,1]| (h=0.001)$ ', 'Interpreter', 'latex');


## exponensial function
import numpy as np
def expfunc(x,a,b):
if x>=a and x<=b:
return np.exp(-x-1)
else:
return 0

# the truncated function f(x) in [-1,1]

def func(x):
if abs(x)<=1:
return 1/(1+x^2)
else:
return 0

# the truncated function f(x) in [a,b]

def linearfunc(x,a,b):

if x>=a and x<=b:

return x
else:
return 0

# the truncated function f(x) in [a,b]

def quadfunc(x,a,b):

if x>=a and x<=b:

return x**2
else:
return 0

# the truncated function f(x) in [a,b]

def rectfunc(x,a,b):

if x>=a and x<=b:

return 1
else:
return 0
# the sinus function f(x) in [a,b]

def sinfunc(x,a,b):

if x>=a and x<=b:

y=np.sin((2*np.pi*(x-a))/(b-a))
else:
return 0
####
def Kw(w,h):
return (np.sin(np.pi*w.*h)./(np.pi*w.*h)).**4*(6/(2*np.cos(2*np.pi*w.*h)+4))

Matlab
% Computation of the sum of the optimal coefficients in the case m=1
%

clear
close all

%
% The interval [a1,b1] is the interval, where functions are given
%
a1=0;
b1=1;
%
% The interval [a,b] is the wider interval than [a1,b1] in which we extent
% the given function with zero value
%
a=0;
b=1;
%
% the set of nodes in the interval [a,b] with respect to variable x and step h=0.1
%
x=a:0.01:b;
%
% the set of nodes in the interval [a1,b1] with respect to variable x and
% step h=0.1
%
%
% defining the step h in the interval [a1,b1] with respect to x
%
%h1=(b1-a1)/N1;
%
% the set of nodes in the interval [a,b] with respect to variable w, with
% step tau=0.01
%
w=a:0.01:b;
M=length(w)-1;
%
% defining the step tau in the interval [a,b] with respect to w
%
h=(b-a)/M;
N=length(x)-1;
%
% f0(n) is the set of values of the rec[a1,b1] function
% f1(n) is the set of values of the sin[a1,b1] function
% f2(n) is the set of values of the exp[a1,b1] function
%

%%
% the function g0, g1 and g2 are Fourier transforms of the functions

% rect[-1,1], sin[-1,1] and exp[-1,1]
% C is coefficients matrix
%
% h=0.01
C=zeros(length(w),length(x));
E0=zeros(1,length(w));
E1=zeros(1,length(w));
E2=zeros(1,length(w));
E11=zeros(1,length(w));
E21=zeros(1,length(w));
f0=zeros(1,length(x));
f1=zeros(1,length(x));
f2=zeros(1,length(x));
f11=zeros(1,length(x));
f21=zeros(1,length(x));
for n=1:(N+1)
f0(n)=rectfunc(x(n),a1,b1);
f1(n)=sinfunc(x(n),a1,b1); %sin
f2(n)=expfunc(x(n),a1,b1); %exp
f11(n)=linearfunc(x(n),a1,b1); %linear
f21(n)=quadfunc(x(n),a1,b1); %quad
end

for m=1:(M+1)

E0(m)=g0(w(m),a1,b1);
E1(m)=gsin(w(m),a1,b1);
E2(m)=gexp(w(m),a1,b1);
E11(m)=g1(w(m),a1,b1);
E21(m)=g2(w(m),a1,b1);
C(m,1)=C0(w(m),h,a);
C(m,N+1)=CN(w(m),h,b);
for n=2:N
C(m,n)=Cb(n-1,w(m),h,a);
end
end
F0=C*f0';
F1=C*f1';
F2=C*f2';
F11=C*f11';
F21=C*f21';
%%
% F0, F1 and F2 are Fourier transforms of rect, linear and quadratic
% functions
%
%
%% h=0.1
w1=a:0.1:b;
M_1=length(w1)-1;
h_1=(b-a)/M_1;
x1=a:0.1:b;
N1=length(x1)-1;

C_1=zeros(length(w1),length(x1));

E0_1=zeros(1,length(w1));
E1_1=zeros(1,length(w1));
E2_1=zeros(1,length(w1));
E11_1=zeros(1,length(w1));
E21_1=zeros(1,length(w1));

f0_1=zeros(1,length(x1));

f1_1=zeros(1,length(x1));
f2_1=zeros(1,length(x1));
f11_1=zeros(1,length(x1));
f21_1=zeros(1,length(x1));
for n=1:(N1+1)
f0_1(n)=rectfunc(x1(n),a1,b1);
f1_1(n)=sinfunc(x1(n),a1,b1); %sin
f2_1(n)=expfunc(x1(n),a1,b1); %exp
f11_1(n)=linearfunc(x1(n),a1,b1); %linear
f21_1(n)=quadfunc(x1(n),a1,b1); %quad
end

for m=1:(M_1+1)

E0_1(m)=g0(w1(m),a1,b1);
E1_1(m)=gsin(w1(m),a1,b1);
E2_1(m)=gexp(w1(m),a1,b1);
E11_1(m)=g1(w1(m),a1,b1);
E21_1(m)=g2(w1(m),a1,b1);
C_1(m,1)=C0(w1(m),h_1,a);
C_1(m,M_1+1)=CN(w1(m),h_1,b);
for n=2:N1
C_1(m,n)=Cb(n-1,w1(m),h_1,a);
end
end
F0_1=C_1*f0_1';
F1_1=C_1*f1_1';
F2_1=C_1*f2_1';
F11_1=C_1*f11_1';
F21_1=C_1*f21_1';

%% h=0.001

w2=a:0.001:b;
M_2=length(w2)-1;
h_2=(b-a)/M_2;
x2=a:0.001:b;
N2=length(x2)-1;

C_2=zeros(length(w2),length(x2));

E0_2=zeros(1,length(w2));
E1_2=zeros(1,length(w2));
E2_2=zeros(1,length(w2));
E11_2=zeros(1,length(w2));
E21_2=zeros(1,length(w2));

f0_2=zeros(1,length(x2));

f1_2=zeros(1,length(x2));
f2_2=zeros(1,length(x2));
f11_2=zeros(1,length(x2));
f21_2=zeros(1,length(x2));

for n=1:(N2+1)

f0_2(n)=rectfunc(x2(n),a1,b1);
f1_2(n)=sinfunc(x2(n),a1,b1); %sin
f2_2(n)=expfunc(x2(n),a1,b1); %exp
f11_2(n)=linearfunc(x2(n),a1,b1); %linear
f21_2(n)=quadfunc(x2(n),a1,b1); %quad
end

for m=1:(M_2+1)

E0_2(m)=g0(w2(m),a1,b1);
E1_2(m)=gsin(w2(m),a1,b1);
E2_2(m)=gexp(w2(m),a1,b1);
E11_2(m)=g1(w2(m),a1,b1);
E21_2(m)=g2(w2(m),a1,b1);
C_2(m,1)=C0(w2(m),h_2,a);
C_2(m,M_2+1)=CN(w2(m),h_2,b);
for n=2:N2
C_2(m,n)=Cb(n-1,w2(m),h_2,a);
end
end
F0_2=C_2*f0_2';
F1_2=C_2*f1_2';
F2_2=C_2*f2_2';
F11_2=C_2*f11_2';
F21_2=C_2*f21_2';

%% []
xx=-2:0.001:2;

ff=zeros(1,length(xx));
fexp=zeros(1,length(xx));
fsin=zeros(1,length(xx));
flinear=zeros(1,length(xx));
fquad=zeros(1,length(xx));
for n=1:length(xx)
ff(n)=rectfunc(xx(n),a1,b1);
fsin(n)=sinfunc(xx(n),a1,b1);
fexp(n)=expfunc(xx(n),a1,b1);
flinear(n)=linearfunc(xx(n),a1,b1);
fquad(n)=quadfunc(xx(n),a1,b1);
end
%% OQF and FFT [f0]
figure %f11
subplot(3,2,[1,2])
plot(xx,ff)
axis([-1.5 1.5 -0.5 1.5])
title('$f_0$','Interpreter','latex');

subplot(3,2,3)

plot(w,real(E0),w1,real(F0_1),w,real(F0'),w2,real(F0_2'))
legend('FT','OQF (h=0.1)','OQF (h=0.01)','OQF (h=0.001)')
axis([a b -5 5])
title('$F[f_0]$','Interpreter','latex');
subplot(3,2,4)
plot(w1,abs(real(F0_1')-real(E0_1)))
title('$|\Re R_{f_{0,\omega}}[-1,1]| (h=0.1)$ ','Interpreter','latex');

subplot(3,2,5)

plot(w,abs(real(F0')-real(E0)))
title('$|\Re R_{f_{0,\omega}}[-1,1]| (h=0.01)$ ','Interpreter','latex');

subplot(3,2,6)

plot(w2,abs(real(F0_2')-real(E0_2)))
title('$|\Re R_{f_{0,\omega}}[-1,1]| (h=0.001)$ ','Interpreter','latex');

%% OQF and FFT [fsin]

figure
subplot(3,2,[1,2])
plot(xx,fsin)
axis([-1.5 1.5 -1.5 1.5])
title('$f_{\sin}$','Interpreter','latex');

subplot(3,2,3)

plot(w,real(E1),w1,real(F1_1),w,real(F1'),w2,real(F1_2'))
legend('FT','OQF (h=0.1)','OQF (h=0.01)','OQF (h=0.001)')
axis([a b -5 5])
title('$F[f_{\sin}]$','Interpreter','latex');
subplot(3,2,4)
plot(w1,abs(real(F1_1')-real(E1_1)))
title('$|\Re R_{f_{{\sin},\omega}}[-1,1]| (h=0.1)$ ','Interpreter','latex');

subplot(3,2,5)

plot(w,abs(real(F1')-real(E1)))
title('$|\Re R_{f_{{\sin},\omega}}[-1,1]| (h=0.01)$ ','Interpreter','latex');

subplot(3,2,6)

plot(w2,abs(real(F1_2')-real(E1_2)))
title('$|\Re R_{f_{{\sin},\omega}}[-1,1]| (h=0.001)$ ','Interpreter','latex');
%% OQF and FFT [fexp]
figure
subplot(3,2,[1,2])
plot(xx,fexp)
axis([-1.5 1.5 -1.5 1.5])
title('$f_{\exp}$','Interpreter','latex');

subplot(3,2,3)

plot(w,real(E2),w1,real(F2_1),w,real(F2'),w2,real(F2_2'))
legend('FT','OQF (h=0.1)','OQF (h=0.01)','OQF (h=0.001)')
axis([a b -5 5])
title('$F[f_{\exp}]$','Interpreter','latex');
subplot(3,2,4)
plot(w1,abs(real(F2_1')-real(E2_1)))
title('$|\Re R_{f_{{\exp},\omega}}[-1,1]| (h=0.1)$ ','Interpreter','latex');

subplot(3,2,5)

plot(w,abs(real(F2')-real(E2)))
title('$|\Re R_{f_{{\exp},\omega}}[-1,1]| (h=0.01)$ ','Interpreter','latex');

subplot(3,2,6)

plot(w2,abs(real(F2_2')-real(E2_2)))
title('$|\Re R_{f_{{\exp},\omega}}[-1,1]| (h=0.001)$ ','Interpreter','latex');

%% OQF and FFT [flinear]

figure
subplot(3,2,[1,2])
plot(xx,flinear)
axis([-1.5 1.5 -1.5 1.5])
title('$f_{linear}$','Interpreter','latex');

subplot(3,2,3)

plot(w,real(E11),w1,real(F11_1),w,real(F11'),w2,real(F11_2'))
legend('FT','OQF (h=0.1)','OQF (h=0.01)','OQF (h=0.001)')
axis([a b -5 5])
title('$F[f_linear]$','Interpreter','latex');
subplot(3,2,4)
plot(w1,abs(real(F11_1')-real(E11_1)))
title('$|\Re R_{f_linear},\omega}}[-1,1]| (h=0.1)$ ','Interpreter','latex');

subplot(3,2,5)

plot(w,abs(real(F11')-real(E11)))
title('$|\Re R_{f_{linear,\omega}}[-1,1]| (h=0.01)$ ','Interpreter','latex');

subplot(3,2,6)

plot(w2,abs(real(F11_2')-real(E11_2)))
title('$|\Re R_{f_linear},\omega}}[-1,1]| (h=0.001)$ ','Interpreter','latex');

%% OQF and FFT [fquad]

figure
subplot(3,2,[1,2])
plot(xx,fquad)
axis([-1.5 1.5 -1.5 1.5])
title('$f_{\quad}$','Interpreter','latex');

subplot(3,2,3)

plot(w,real(E21),w1,real(F21_1),w,real(F21'),w2,real(F21_2'))
legend('FT','OQF (h=0.1)','OQF (h=0.01)','OQF (h=0.001)')
axis([a b -5 5])
title('$F[f_{\quad}]$','Interpreter','latex');
subplot(3,2,4)
plot(w1,abs(real(F21_1')-real(E21_1)))
title('$|\Re R_{f_{{\quad},\omega}}[-1,1]| (h=0.1)$ ','Interpreter','latex');

subplot(3,2,5)

plot(w,abs(real(F21')-real(E21)))
title('$|\Re R_{f_{{\exp},\omega}}[-1,1]| (h=0.01)$ ','Interpreter','latex');

subplot(3,2,6)

plot(w2,abs(real(F21_2')-real(E21_2)))
title('$|\Re R_{f_{{\quad},\omega}}[-1,1]| (h=0.001)$ ','Interpreter','latex');
%
% C_0 is the first coefficient of the optimal quadrature formula
%
%
function y=C0(w,h,a)
if w==0
y=h/2;
else
y=(h/2)*Kw(w,h).*exp(2*pi*i*w.*a);
end
end

%
% C_0 is the first coefficient of the optimal quadrature formula

%
%
function y=Cb(beta,w,h,a)
if w==0
y=h;
else
y=h*Kw(w,h).*exp(2*pi*i*w.*(h*beta+a));
end
end

%
% C_N is the last coefficient of the optimal quadrature formula

%
%
function y=CN(w,h,b)
if w==0
y=h/2;
else
y=(h/2)*Kw(w,h).*exp(2*pi*i*w.*b);
end
end

function y=expfunc(x,a,b)

if x>=a && x<=b

y=exp(-x-1);
else
y=0;

end

%
%
% the truncated function f(x) in [-1,1]
%
function y=func(x)

if abs(x)<=1

y=1/(1+x^2);
else
y=0;
end

%
%
% the truncated function f(x) in [a,b]

%
function y=linearfunc(x,a,b)

if x>=a && x<=b

y=x;
else
y=0;
end
%
%
% the truncated function f(x) in [a,b]
%
function y=quadfunc(x,a,b)

if x>=a && x<=b

y=x^2;
else
y=0;
end

%
%
% the truncated function f(x) in [a,b]

%
function y=rectfunc(x,a,b)

if x>=a && x<=b

y=1;
else
y=0;
end
end

% the sinus function f(x) in [a,b]

%
function y=sinfunc(x,a,b)

if x>=a && x<=b

y=sin((2*pi*(x-a))/(b-a));
else
y=0;

end

%
% Fourier transform of f(x)=1 denoted by g0(w)
%
% g0(w,a,b)=\int_a^b 1 [exp(2\pi iwx)] dx
%

function y=g0(w,a,b)

if w==0
y=b-a;

else
y=-i*(exp(2*pi*i*w*b)-exp(2*pi*i*w*a))/(2*pi*w);
end
end

%
% Fourier transform of f(x)=x denoted by g1(w)

%
% g1(w,a,b)=\int_a^b [x exp(2\pi iwx)]dx
%

function y=g1(w,a,b)

if w==0
y=(b^2-a^2)/2;

else
y=(exp(2*pi*i*w*b)*b-exp(2*pi*i*w*a)*a)/(2*pi*i*w)-(exp(2*pi*i*w*b)-exp(2*pi*i*w*a))/((2*pi*i*w)^2);
end
%
% Fourier transform of f(x)=x^2 denoted by g2(w)
%
% g2(w,a,b)=\int_a^b [x^2 exp(2\pi iwx)]dx
%

function y=g2(w,a,b)

if w==0
y=(b^3-a^3)/3;

else
y=(exp(2*pi*i*w*b)*b^2-exp(2*pi*i*w*a)*a^2)/(2*pi*i*w)...
-2*(exp(2*pi*i*w*b)*b-exp(2*pi*i*w*a)*a)/((2*pi*i*w)^2)...
+2*(exp(2*pi*i*w*b)-exp(2*pi*i*w*a))/((2*pi*i*w)^3);
End

function y=gexp(w,a,b)

% y=(a-b)*(exp(-1+2*i*pi*b*w)-exp(2*i*pi*a*w))*(1-2*i*pi*a*w+2*i*pi*b*w)/...
% (4*pi^2*a^2*w^2-8*pi^2*a*b*w^2+4*pi^2*b^2*w^2+1);
%y=exp(2*pi*i*b*w-b)/(2*i*pi*w-1)-exp(2*pi*i*a*w-a)/(2*i*pi*w-1);
y=exp(2*pi*i*b*w-b-1)/(2*i*pi*w-1)-exp(2*pi*i*a*w-a-1)/(2*i*pi*w-1);

end
function y=gsin(w,a,b)

y=(a-b)*(exp(2*1i*pi*a*w)-exp(2*1i*pi*b*w))/(2*pi*(a^2*w^2-2*a*b*w^2+b^2*w^2-1));
%y=(exp(2*i*pi*a*w)*(cos(pi*a)-2*i*w*sin(pi*a))+exp(2*i*pi*b*w)*(-cos(pi*b)+2*i*w*sin(pi*b)))/(pi*(4*w^2-1));
% if a=-b then
% y=2*i*a*sin(2*pi*a*w)/(pi*(4*a^2*w^2-1))

end

function y=Kw(w,h)
y=(sin(pi*w.*h)./(pi*w.*h)).^4*(6/(2*cos(2*pi*w.*h)+4));
end