Homec4science
Diffusion GYACOMO (OBSOLETE: LATEST VERSION ON gitlab.epfl.ch) (divided_output) Browse

gyacomo/matlab/plot/plot_kernels.m
076e4b173ededivided_output

Last Change
Raw File

plot_kernels.m
View Options

if 0
%% Kernels
kmax=5;
nmax=6;
kx = linspace(0,kmax,100);
figure
for n_ = 0:nmax
plot(kx,kernel(n_,kx),'DisplayName',['$\mathcal{K}_{',num2str(n_),'}$']);hold on;
end
ylim_ = ylim;
plot(kx(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
% J0 = besselj(0,kx);
% plot(kx,J0,'-r','DisplayName','$J_0$');
legend('show'); xlabel('k'); title('Kernel functions $\mathcal{K}_n$');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 0
%% Bessels and approx
vperp = linspace(0,1.5,4);
nmax1=5;
nmax2=10;
kmax=7;
figure
for i = 1:4
subplot(2,2,i)
v_ = vperp(i);
kx = linspace(0,kmax,100);
J0 = besselj(0,kx*v_);
A1 = 1 - kx.^2*v_^2/4;
K1 = zeros(size(kx));
K2 = zeros(size(kx));
for n_ = 0:nmax1
K1 = K1 + kernel(n_,kx).*polyval(LaguerrePoly(n_),v_^2);
end
for n_ = 0:nmax2
K2 = K2 + kernel(n_,kx).*polyval(LaguerrePoly(n_),v_^2);
end
plot(kx,J0,'-k','DisplayName','$J_0(kv)$'); hold on;
plot(kx,A1,'-r','DisplayName','$1 - k^2 v^2/4$');
plot(kx,K1,'--b','DisplayName',['$\sum_{n=0}^{',num2str(nmax1),'}\mathcal K_n(k) L^n(v)$']);
plot(kx,K2,'-b','DisplayName',['$\sum_{n=0}^{',num2str(nmax2),'}\mathcal K_n(k) L^n(v)$']);
ylim_ = [-0.5, 1.0];
plot(kx(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
ylim(ylim_); xlabel('$k$')
legend('show'); grid on; title(['$v = ',num2str(v_),'$'])
end
%%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 0
%% from Sum_n Kn x Ln x Lj to Sum_n Kn Sum_s dnjs x Ls
vperp = [2];
kx = linspace(0,10,200);
Jmax = 15;
j = 10;
fig = figure; set(gcf, 'Position', [100, 100, numel(vperp)*300, 250])
% suptitle(['$J_{max}=',num2str(Jmax),', j=',num2str(j),'$'])
Kn = @(x__,y__) kernel(x__,y__);
for i = 1:numel(vperp)
subplot(1,numel(vperp),i)
v_ = vperp(i);
% J0 x Lj
J0Lj = besselj(0,kx*v_)*polyval(LaguerrePoly(j),v_^2);
% Taylor exp of J0
A1 = (1 - kx.^2*v_^2/4)*polyval(LaguerrePoly(j),v_^2);
% Sum_n^Nmax Kn x Ln x Lj
KLL = zeros(size(kx));
for n_ = 0:Jmax
KLL = KLL + Kn(n_,kx).*polyval(LaguerrePoly(n_),v_^2);
end
KLL = KLL.*polyval(LaguerrePoly(j),v_^2);
% Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
KdL1 = zeros(size(kx));
for n_ = 0:Jmax-j
tmp_ = 0;
for s_ = 0:n_+j
tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2);
end
KdL1 = KdL1 + Kn(n_,kx).*tmp_;
end
% Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
KdL2 = zeros(size(kx));
for n_ = 0:Jmax
tmp_ = 0;
for s_ = 0:min(Jmax,n_+j)
tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2);
end
KdL2 = KdL2 + Kn(n_,kx).*tmp_;
end
% plots
plot(kx,J0Lj,'-k','DisplayName','$J_0 L_j$'); hold on;
plot(kx,KLL,'-','DisplayName',['$\sum_{n=0}^{J}\mathcal K_n L^n L_j$']);
plot(kx,KdL1,'-.','DisplayName',['$\sum_{n=0}^{J-j}\mathcal K_n \sum_{s=0}^{n+j} d_{njs} L^s$']);
ylim_ = ylim;
plot(kx,A1,'-','DisplayName','$(1 - k^2 v^2/4) L_j$'); hold on;
plot(kx,KdL2,'--','DisplayName',['$\sum_{n=0}^{J}\mathcal K_n \sum_{s=0}^{\min(J,n+j)} d_{njs} L^s$']);
% plot(kx(end)*[2/3 2/3],ylim_,'--k','DisplayName','AA');
ylim(ylim_);
xlabel('$k_{\perp}$')
legend('show');
grid on; title(['$v = ',num2str(v_),'$',' $J_{max}=',num2str(Jmax),', j=',num2str(j),'$'])
end
FIGNAME = ['/home/ahoffman/Dropbox/Applications/Overleaf/Hermite-Laguerre Z-pinch (HeLaZ Doc.)/figures/J0Lj_approx_J','_',num2str(Jmax),'_j_',num2str(j),'.eps'];
% saveas(fig,FIGNAME,'epsc');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 1
%% Convergence analysis
kx = linspace(0,10,200);
v_A = linspace(0,1,10);
J_A = 1:15;
[YY,XX] = meshgrid(v_A,J_A);
klimLL = zeros(size(XX));
klimdL1= zeros(size(XX));
Kn = @(x__,y__) kernel(x__,y__);
fig = figure; set(gcf, 'Position', [100, 100, 500, 400]);
for j = 5
for ij_ = 1:numel(J_A)
iv_ = 1;
for v_ = v_A
eps = abs(0.01*polyval(LaguerrePoly(j),v_^2));
Jmax = J_A(ij_);
% J0 x Lj
J0Lj = besselj(0,kx*v_)*polyval(LaguerrePoly(j),v_^2);
% Taylor exp of J0
A1 = (1 - kx.^2*v_^2/4)*polyval(LaguerrePoly(j),v_^2);
% Sum_n^Nmax Kn x Ln x Lj
KLL = zeros(size(kx));
for n_ = 0:Jmax
KLL = KLL + Kn(n_,kx).*polyval(LaguerrePoly(n_),v_^2);
end
KLL = KLL.*polyval(LaguerrePoly(j),v_^2);
% Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
KdL1 = zeros(size(kx));
for n_ = 0:Jmax-j
tmp_ = 0;
for s_ = 0:n_+j
tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2);
end
KdL1 = KdL1 + Kn(n_,kx).*tmp_;
end
% Sum_n^Nmax Kn Sum_s^Smax dnjs x Ls
KdL2 = zeros(size(kx));
for n_ = 0:Jmax
tmp_ = 0;
for s_ = 0:min(Jmax,n_+j)
tmp_ = tmp_ + dnjs(n_,j,s_)*polyval(LaguerrePoly(s_),v_^2);
end
KdL2 = KdL2 + Kn(n_,kx).*tmp_;
end
% errors
err_ = abs(J0Lj-KLL);
idx_ = 1;
while(err_(idx_)< eps && idx_ < numel(err_))
idx_ = idx_ + 1;
end
klimLL(ij_,iv_) = kx(idx_);
err_ = abs(J0Lj-KdL1);
idx_ = 1;
while(err_(idx_)< eps && idx_ < numel(err_))
idx_ = idx_ + 1;
end
klimdL1(ij_,iv_) = kx(idx_);
iv_ = iv_ + 1;
end
end
% plot
% plot(JmaxA,klimLL,'o-'); hold on
% plot(J_A,klimdL1,'o-'); hold on
end
surf(XX,YY,klimdL1)
xlabel('$J_{max}$'); ylabel('$v_{\perp}$');
title(['$k$ s.t. $\epsilon=1\%$, $j=',num2str(j),'$'])
ylim([0,1]); grid on;
end
if 0
%% Test Lj Ln = sum_s dnjs Ls
j = 10;
n = 10;
smax = n+j;
vperp = linspace(0,5,20);
LjLn = zeros(size(vperp));
dnjsLs = zeros(size(vperp));
dnjsLs_bin = zeros(size(vperp));
dnjsLs_stir = zeros(size(vperp));
i=1;
for x_ = vperp
LjLn(i) = polyval(LaguerrePoly(j),x_)*polyval(LaguerrePoly(n),x_);
dnjsLs(i) = 0;
dnjsLs_bin(i) = 0;
dnjsLs_stir(i) = 0;
for s_ = 0:smax
dnjsLs(i) = dnjsLs(i) + dnjs(n,j,s_)*polyval(LaguerrePoly(s_),x_);
dnjsLs_bin(i) = dnjsLs_bin(i) + dnjs_fact(n,j,s_)*polyval(LaguerrePoly(s_),x_);
dnjsLs_stir(i) = dnjsLs_stir(i) + dnjs_stir(n,j,s_)*polyval(LaguerrePoly(s_),x_);
end
i = i+1;
end
figure
plot(vperp,LjLn); hold on;
plot(vperp,dnjsLs,'d')
plot(vperp,dnjsLs_bin,'o')
plot(vperp,dnjsLs_stir/dnjsLs_stir(1),'x')
% We see that starting arround j = 18, n = 0, the stirling formula is the only
% approximation that does not diverge from the target function. However it
% performs badly for non 0 n...
end
c4science · Help