%%% Linear Computational Instability for Diffusion Equation % % d_phi/dt = kappa d^2 phi/dx^2 (1) % % The general finite difference equation is % % (phi_j - phi_j^0)/dt = kappa { (2) % theta [ phi_(j+1) - 2 phi_j + phi_j-1) ]/ dx^2 % + (1-alpha) [ phi_(j+1)^0 - 2 phi_j^0 + phi_j-1)^0 ]/ dx^2 % } % where phi_j are at future time, phi_j^0 are known values % % Take spatial Discrete Fourier Transform (assume infinite domain) % % PHI(k) = sum_{j=-inf}^inf [ exp(ik jdx) phi_j ] dx % % Take DFT of (2) and shift dummy indices j+1 -> m, j-1-> n % (recall that the sum is from -inf to inf) to recover standard DFT % from j+1 and j-1 terms. % % PHI(k) - PHI^0(k) = (kappa dt/dx^2) { % theta [ PHI(k) exp(-ikdx) - 2 PHI(k) + PHI(k) exp(ikdx) ] % + (1-alpha)[ PHI_0(k) exp(-ikdx) - 2 PHI^0(k) + PHI^0(k) exp(ikdx)] % } % Collecting terms % % PHI(k) = PHI^0(k) F(k) % % where F(k) is the transfer function for a time step dt % for wavenumber k % % F(k) = [ 1 - 2(1-alpha)(kappa dt/dx^2) (1 - cos(k dx) ) ] / ... % [ 1 + 2 alpha (kappa dt/dx^2) (1 - cos(k dx) ) ] % % The Nyquist wavenumber k_N = pi/dx % is highest wavenumber representable on a uniform grid % ------------------------------------------------------- % % Also demonstrate linear computational instability in time domain % % Ed Waddington 2008-05-11 % -------------------------------------------------------------------- % % set up diffusivity k kappa = 1; % % set up grid length Lx Lx = 20; % % set up discretization dx, and nodes vector x_vec dx = 1.0; x_vec = 0:dx:Lx; nx = length(x_vec); % % % set duration Lt of time stepping % Lt = 10; % % set time step dt = input('OK, what is dt for this test? '); % nt = fix( Lt/dt ) + 1; % % print out stability parameter beta = kappa dt/dx^2 % ---------------------------------------------------- beta = 2*kappa*dt/dx^2; stability_string = num2str( 2*kappa*dt/dx^2); % message_string ='beta = [ 2 kappa dt/dx^2 ] = '; panic_message = [ message_string stability_string ]; % disp( panic_message ) % % Set time-step mixing fraction alpha for spatial terms % i.e. alpha*phi + (1-alpha)*phi_0 % alpha = 1 Fully Implicit % = 0.5 Crank-Nicolson % = 0 Explicit % ------------------------------ % % alpha should exceed (beta-1)/2 beta % ---------------------------------------- alpha_test_string = ('For stability: choose alpha > (beta - 1) / (2 beta) = '); alpha_test_message = [ alpha_test_string num2str( (beta - 1)/(2*beta) ) ]; disp( alpha_test_message ); % alpha = input('OK what value for time-partition parameter alpha? '); % alpha_label_string ='alpha = '; alpha_value_string =num2str( alpha ); alpha_message = [ alpha_label_string alpha_value_string ]; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Now set up wavenumber parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % k_N = pi/dx; % Nyquist wavenumber % n_k = 51; % number of points in wavenumber domain % dk = k_N/(n_k - 1); % k_vec = 0: dk: k_N; % % Evaluate transfer function for continuous diffusion equation % ------------------------------------------------------------- F_analytic = exp(- kappa * k_vec.^2 * dt ); % % Evaluate transfer function for discrete diffusion equation % ------------------------------------------------------------ F = ( 1 - (1-alpha)*(2*kappa*dt/dx^2)*(1 - cos(k_vec*dx) ) ) ./ ... ( 1 + alpha*(2*kappa*dt/dx^2)*(1 - cos(k_vec*dx) ) ); % % show kappa*dt/dx^2 % ------------------- % disp('2*kappa*dt/dx^2 = '); disp(2*kappa*dt/dx^2) % % Plot transfer functions F_analytical(k) and F(k) % ------------------------------------------------- % figure(1) plot(k_vec, F_analytic, 'Color', 'k' ) hold on;grid on;box on plot( k_vec, zeros(size(k_vec)), 'linewidth', 2 ) xlabel('Wavenumber k') ylabel('Transfer function F(k) for time-step dt') % % put white patch to cover old dt and alpha messages patch( [ 0.05*pi 0.05*pi 0.2*pi 0.2*pi ], ... [ 0.1 0.25 0.25 0.1 ], 'w' ) % % write current dt and alpha onto plot dt_string = ['dt = ' num2str(dt) ]; text(0.06*pi,0.22, dt_string ) % alpha_string = ['\alpha = ' num2str(alpha) ]; text(0.06*pi,0.12, alpha_string ) % title('Black - Analytical Transfer function for time step dt') % color = input('What color of line to plot _this_ time? '); plot(k_vec, F, 'Color', color ) % % % Plot transfer function for repeated steps to give % total elapsed time t=1 % When F(k) < 0, and 1/dt is not an integer, % the product of transfer functions needed to represent passage % of t=1 may be complex. However, I plot only the magnitude. % --------------------------------------------------------------- figure(2) plot(k_vec, F_analytic.^(1/dt), 'Color', 'k' ) hold on;grid on;box on xlabel('Wavenumber k') title('Multi-step Transfer function F(k) to simulate t=1') ylabel('Transfer function F(k)') text(0.7*pi,0.85, 'repeated to simulate t=1' ) % plot(k_vec, abs( F.^(1/dt) ), 'Color', color ) % % % put white patch to cover old dt and alpha messages patch( [ 0.02*pi 0.02*pi 0.2*pi 0.2*pi ], ... [ 0.1 0.4 0.4 0.1 ], 'w' ) % % write current dt and alpha onto plot dt_string = ['dt = ' num2str(dt) ]; text(0.06*pi,0.32, dt_string ) % alpha_string = ['\alpha = ' num2str(alpha) ]; text(0.06*pi,0.12, alpha_string ) % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % end of wavenumber calculation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Now find finite-difference solution in time domain %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % How long to run simulation? % set duration Lt of time stepping and number nt of time steps % ------------------------------------------------------------ Lt = input('How long to run this simulation? '); nt = fix( Lt/dt ) + 1; % % set string for plot-ends time time_ends = [ ' Run ends at ' num2str(Lt) ]; % % % set up initial condition (sum of a few sinusoids) % The b_N term puts in a little ripple at Nyquist wavenumber. % Does that sinificantly "feed" runaway growth? % --------------------------------------------------- a_0 = 10; a_1= 1.0; b_1 = 2.0; b_2 = 2.0; b_3 = 2.0; b_4 = 2.0; b_N = 0.0; phi_P_0 = a_0 + a_1 * x_vec/Lx ... + b_1 * sin( pi*x_vec/Lx ) ... + b_2 * sin(2*pi*x_vec/Lx) ... + b_3 * sin(3*pi*x_vec/Lx) ... + b_4 * sin(4*pi*x_vec/Lx) ... + b_N * cos((nx-1)*pi*x_vec/Lx); % % % Boundary Conditions % -------------------- % Set to hold steady from initial condition. % % specified value at x = 0 phi_0 = a_0 + b_N; % % specified gradient at x = Lx phi_prime_Lx = a_1/Lx ... + b_1 * (pi/Lx) * cos( pi*x_vec(end)/Lx ) ... + b_2 * (2*pi/Lx) * cos(2*pi*x_vec(end)/Lx ) ... + b_3 * (3*pi/Lx) * cos(3*pi*x_vec(end)/Lx ) ... + b_4 * (4*pi/Lx) * cos(4*pi*x_vec(end)/Lx ) ... + b_N * (nx*pi/Lx) * cos(nx*pi*x_vec(end)/Lx ); % % % set up shifted phi arrays for phi_E and phi_W phi_D_0 = [ phi_P_0(2:end) 0 ]; phi_U_0 = [ 0 phi_P_0(1:end-1) ]; % % % set up coefficient vectors a_D = repmat( (alpha*kappa*dt/dx^2), nx, 1)'; a_U = a_D; a_P = ( 1 + a_D + a_U ); % a_E = zeros( size(a_P) ); a_W = a_E; % a_D_0 = repmat( ( (1-alpha)*kappa*dt/dx^2), nx, 1)'; a_U_0 = a_D_0; a_P_0 = 1 - a_D_0 - a_U_0; % b = a_P_0 .*phi_P_0 + a_D_0 .*phi_D_0 + a_U_0 .*phi_U_0; % % % Introduce Boundary Conditions into equations % % fixed value at x_vec(1) a_P(1) = 1; a_U(1) = 0; a_D(1) = 0; % a_P_0(1) = 0; a_U_0(1) = 0; a_D_0(1) = 0; b(1) = phi_0; % % fixed gradient at x_vec(end) a_P(end) = 1/dx; a_U(end) = -1/dx; a_D(end) = 0; % a_P_0(end) = 0; a_U_0(end) = 0; a_D_0(end) = 0; b(end) = phi_prime_Lx; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % end of setup for time domain %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %S % % Start plot with initial-condition profile fig_no = input('Which figure (>=3) to use for finite-difference solution? '); figure(fig_no) % % set y scale limits on plot phi_max = 25; % plot(x_vec, phi_P_0, 'k') hold on;grid on;box on; xlabel('Depth z') ylabel('Solution \phi(z)') text( 0.1*Lx, 0.95*phi_max, panic_message ) text( 0.1*Lx, 0.90*phi_max, alpha_message ) text( 0.1*Lx, 0.85*phi_max, alpha_test_message); text( 0.1*Lx, 0.80*phi_max, time_ends ) title('Stability tests with Finite Differences') set(gca,'ylim', [ 0, phi_max] ) % pause % % Here it is! Big time-step loop! % --------------------------------- for j = 1:nt-1 % % First, update RHS with current "old" solution % set up shifted phi arrays for phi_0_E and phi_0_W phi_D_0 = [ phi_P_0(2:end) 0 ]; phi_U_0 = [ 0 phi_P_0(1:end-1) ]; % % Update RHS vector, which depends on phi_0 b = a_P_0 .*phi_P_0 + a_D_0 .*phi_D_0 + a_U_0 .*phi_U_0; % % Re-introduce Boundary conditions into RHS vector b(1) = phi_0; b(end) = phi_prime_Lx; % % a_P(end) = 1/dx; a_U(end) = 1/dx; a_D(end) = 0; % a_P_0(end) = 0; a_U_0(end) = 0; a_D_0(end) = 0; % % Lots of equations? Let's solve the suckers all together! phi = solver( a_E, a_W, a_U, a_D, a_P, b ); % % update "old" time step phi_P_0 = phi; % % Add current time step to plot figure(fig_no) plot(x_vec, phi, 'Color', color ) % % show time step and elapsed time time_message = ['Current Time = ' num2str( j*dt) ]; patch( Lx*[0.1 0.1 0.4 0.4], phi_max*[0.73 0.77 0.77 0.73 ], 'w' ) text(0.1*Lx, 0.75*phi_max, time_message) % % pause % end % % % print maximum value of phi - is it going wild? disp('Maximum of |phi| = '); disp( max( abs(phi) ) ); % % %