Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/matlab/fourier_smooth.m
f1da25c60953feature/refactoring

Last Change
Raw File

fourier_smooth.m
View Options

%
% @file fourier_smooth.m
%
% @brief
%
% @copyright
% Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
% SPC (Swiss Plasma Center)
%
% spclibs is free software: you can redistribute it and/or modify it under
% the terms of the GNU Lesser General Public License as published by the Free
% Software Foundation, either version 3 of the License, or (at your option)
% any later version.
%
% spclibs is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.
%
% @authors
% (in alphabetical order)
% @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
%
% Brute force computation of \mu = MAX |G| for
% damped Jacobi and Gauss-Seidel relaxation
%
% Find optimal damping factor \omega for Jacobi relaxations
%
if ~exist('tau','var'), tau=0; end
if ~exist('alpha','var'), alpha=1.0; end
dth=0.01;
theta1=-pi+dth:dth:pi;
theta2=-pi+dth:dth:pi;
nth1=length(theta1);
nth2=length(theta2);
sint12=sin(theta1./2).^2;
sint22=sin(theta2./2).^2;
sint1=sin(theta1);
sint2=sin(theta2);
[S1,S2]=meshgrid(sint1,sint2);
ctau=0.25*alpha*tau;
omega=0.5:0.002:1;
n=length(omega);
G = zeros(nth1,nth2);
for i=1:n
c=2*omega(i)/(1+alpha^2);
for i1=1:nth1
for i2=1:nth2
if or(abs(theta1(i1))>= pi/2, abs(theta2(i2)) >= pi/2);
G(i1,i2) = abs(1-c*(sint12(i1)+alpha^2*sint22(i2) + ...
ctau*sint1(i1)*sint2(i2)));
end
end
end
[gmax,imax]=max(G);
[mu(i),jmax]=max(gmax);
theta1_opt(i) = theta1(imax(jmax));
theta2_opt(i) = theta2(jmax);
end
[mu_min,i_min]=min(mu);
omega_opt=omega(i_min);
str_title=sprintf(['omega = %.3f, mu = %.3f, alpha = %.2f, tau = ' ...
'%.1f'], omega_opt, mu_min, alpha, tau);
figure
subplot(211)
plot(omega,mu,'LineWidth',2);
xlabel('\omega'); ylabel('\mu')
grid on
title(str_title);
subplot(212)
plot(omega, theta1_opt, omega, theta2_opt,'LineWidth',2);
legend('\theta_{1opt}', '\theta_{2opt}')
xlabel('\omega'); ylabel('optimum \theta')
grid on
%
% \mu for Gauss-Seidel relaxation
%
Ggs = zeros(nth1,nth2);
exp1=complex(cos(theta1),sin(theta1));
exp2=complex(cos(theta2),sin(theta2));
c=2*(1+alpha^2);
ctau=complex(alpha^2, (0.5*alpha*tau).*sint1);
for i1=1:nth1
for i2=1:nth2
if or(abs(theta1(i1))>= pi/2, abs(theta2(i2)) >= pi/2);
Num = exp1(i1) + ctau(i1)*exp2(i2);
Ggs(i1,i2) = abs( Num/(c-conj(Num)) );
end
end
end
[gmax,imax]=max(Ggs);
[mugs,jmax]=max(gmax);
theta1_gs_opt = theta1(imax(jmax));
theta2_gs_opt = theta2(jmax);
fprintf('alpha = %.2f, tau = %.1f, Ggs = %.4f, theta1 = %.4f, theta2 = %.4f\n', alpha, tau, ...
mugs, theta1_gs_opt, theta2_gs_opt);
subplot(211)
hold on
plot(omega,repmat(mugs,1,length(omega)),'r--','LineWidth',2)
c4science · Help