% makeplotS.m % % call plotS to plot the region of absolute stability S % (and the Order Star if plotOS=1 below) for a specific method. % % From http://www.amath.washington.edu/~rjl/fdmbook/chapter7 (2007) method = 'pade12'; % set to desired method, see choices below plotOS = 0; % set to 1 to also plot Order Star disp(' ' ) if plotOS disp(['Plotting S and OS for method ',method]) else disp(['Plotting S for method ',method]) end switch method case 'euler' R = @(z) 1 + z; axisbox = [-3 3 -3 3]; case 'taylor' % Taylor series method of order n n = 12; disp([' order = ' num2str(n)]) R = @(z) taylorexp(n,z); axisbox = [-5 5 -5 5]; if n>7 xy = 0.7*n; axisbox = [-xy xy -xy xy]; end case 'trbdf2' R = @(z) (1 + 5/12 * z) ./ (1 - 7/12 * z + 1/12 * z.^2); axisbox = [-5 15 -10 10]; case 'pade12' R = @(z) (1 + 1/3 * z) ./ (1 - 2/3 * z + 1/6 * z.^2); axisbox = [-2 8 -5 5]; case 'pade22' R = @(z) (1 + 1/2 * z + 1/12 * z.^2) ./ (1 - 1/2 *z + 1/12 * z.^2); axisbox = [-5 5 -5 5]; case 'pade13' R = @(z) (1 + 1/4 * z) ./ (1 - 3/4 * z + 1/4 * z.^2 - 1/24 * z.^3); axisbox = [-3 7 -5 5]; case 'sdirk3' gamma = 2/3; R = @(z) (1 + z*(1-2*gamma) + z.^2 * (1/2 - 2*gamma + gamma^2)) ./ ... (1 - gamma*z).^2; axisbox = [-5 35 -20 20]; case 'radau5' R = @(z) (1 + 2/5 * z + 1/20 * z.^2) ./ ... (1 - 3/5 * z + 3/20 * z.^2 - 1/60 * z.^3); axisbox = [-5 15 -10 10]; otherwise disp(' *** method not recognized') end % Now that R is set, plot S: plotS(R,plotOS,axisbox)