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