function m20121024_01_solveFulopHelanderODE() % 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 = 0.7; gamma = 0.4; % Number of theta gridpoints for linear system N = 15; % 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; % 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; % Plot results: figure(1) clf numRows=2; numCols=1; % Number of point to use for interpolation: NFine = 50; boundary=1; [thetaFine,bFine] = m20120528_03_FourierInterpolation(bs, 0, 2*pi, NFine, boundary); subplot(numRows,numCols,1) plot(thetaFine,bFine) hold on plot(theta,bs,'.r') xlim([0,2*pi]) xlabel('\theta') ylabel('b') [thetaFine,nFine] = m20120528_03_FourierInterpolation(n, 0, 2*pi, NFine, boundary); subplot(numRows,numCols,2) plot(thetaFine,nFine) hold on plot(theta,n,'.r'); xlim([0,2*pi]) xlabel('\theta') ylabel('n') 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