function m20121024_02_solveFulopHelanderODEScan()
% This script solves the top line of equation (9) in
% Helander, Phys Plasmas 5, 3999 (1998)
% in the trace impurity limit alpha=0:
%
% d n/d theta = g * [n - b^2 + gamma * (n - ) * b^2]
%
% The same equation appears as equation (27) in
% Fulop & Helander, Phys Plasmas 8, 3305 (2001).
%
% The script works for an arbitrary normalized magnetic field strength
% B(theta), which is specified as a function at the end of this file.
%
% To solve the linear equation, we first isolate the inhomogeneous term on
% the right-hand side:
%
% n + gamma * (n - ) * b^2 - (1/g) * (d n/d theta) = b^2
%
% We discretize n on a uniformly spaced grid in theta, giving a column
% vector representing n on the grid points. Then notice the left-hand side
% of the above equation is a linear operation applied to n, in particular,
% an integro-differential operation. Since the operator is linear, we may
% represent it as a matrix M. For the derivative term, we use a
% spectral differentiation matrix. The integral term
% is performed by noting that integration is essentially
% multiplication of the n column vector by a row vector. The ODE is solved
% by solving the linear system M * n = b^2 with b^2 represented as a column
% vector.
% Physics parameters:
g_max = 5;
gamma_max = 3;
gs = linspace(0,g_max,17);
gammas = linspace(0,gamma_max,15);
% Number of theta gridpoints for linear system.
N = 16;
% N must be even so there will be a grid point at theta=pi.
if mod(N,2) ~= 0
error('N must be even')
end
halfN = N/2;
% Call an external function to make the Fourier spectral differentiation
% matrix for a uniform grid, as well as the grid itself:
[theta, ddtheta] = fourdif(N,1);
dtheta = theta(2)-theta(1);
% To integrate in theta, sum elements and multiply by dtheta.
Bs = B(theta);
fluxSurfaceAverageB2 = sum(Bs.^2)*dtheta/(2*pi);
bs = Bs ./ sqrt(fluxSurfaceAverageB2);
bSquared = bs.^2;
numGs = numel(gs);
numGammas = numel(gammas);
asymmetry = zeros(numGammas,numGs);
% Loop over all desired values of g and gamma:
for iG = 1:numGs
g = gs(iG);
if g==0
% When g=0, the original ODE is singular, so treat this case
% separately:
asymmetry(:,iG) = 1;
else
for iGamma = 1:numGammas
gamma = gammas(iGamma);
% Now begin construction of the matrix.
% Begin with the identity matrix, representing the n term:
matrix = eye(N);
% Add gamma * n * b^2 term:
matrix = matrix + gamma * diag(bSquared);
% Add -gamma * * b^2 term:
matrix = matrix - gamma * bSquared * (dtheta /(2*pi)* bSquared');
% Add -d/dtheta term:
matrix = matrix - (1/g)*ddtheta;
% We're now finished building the matrix.
% Now the key step: solve the linear system.
n = matrix \ bSquared;
nOutboard = n(1);
nInboard = n(halfN);
asymmetry(iGamma,iG) = nInboard / nOutboard;
end
end
end
% Plot results:
numContours = 15;
figure(1)
clf
contourf(gammas,gs,asymmetry',numContours,'EdgeColor','none')
colorbar
ylabel('g')
xlabel('\gamma')
title('Asymmetry: n(\theta=\pi) / n(theta=0)')
% Compare the numerical calculation to the theoretical limit for g >> 1:
figure(2)
clf
b2In = bSquared(halfN);
b2Out = bSquared(1);
plot(gammas, asymmetry(:,end))
hold on
nOut = b2Out./(1+gammas*b2Out);
nIn = b2In ./(1+gammas*b2In);
plot(gammas, nIn ./ nOut,':r')
legend(['numerical, g=',num2str(gs(end))], 'analytical, g = Inf')
xlabel('\gamma')
ylabel('asymmetry')
title('Benchmark for large g')
end
function BB = B(theta)
% User-supplied function for the field strength |B|.
% You can normalize B any way you like, since
% the output of this function will be divided by **^{-1/2}.
epsilon = 0.3;
BB = 1./(1 + epsilon * cos(theta));
end**