1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 | % Example explicit numerical integration of a linear ODE. % Butcher tableaus (only a and b) for the three methods. tabs = struct(); tabs.euler = struct('a', 0, 'b', 1); tabs.heun = struct('a', [0, 0; 1, 0], 'b', [0.5, 0.5]); tabs.rk4 = struct('a', [0, 0, 0, 0; 0.5, 0, 0, 0; 0, 0.5, 0, 0; 0, 0, 1, 0], ... 'b', [1/6, 1/3, 1/3, 1/6]); tabnames = fieldnames(tabs); % Run each method. A = [0, 1; -1, 0]; f = @(x) A*x; x0 = [1; 0]; T = 2*pi(); Meval = 32; x = struct(); for n = 1:length(tabnames) n = tabnames{n}; Nstep = ceil(Meval/length(tabs.(n).b)); % Number of integration steps. h = T/Nstep; x.(n) = zeros(2, Nstep + 1); x.(n)(:,1) = x0; for i = 1:Nstep x.(n)(:,i + 1) = butcherintegration(x.(n)(:,i), f, h, ... tabs.(n).a, tabs.(n).b); end end % Add exact solution. theta = linspace(0, 2*pi(), 257); x.exact = [cos(theta); -sin(theta)]; tabnames = [{'exact'}; tabnames]; % Make a plot. errorplot2d(x, T, 'SouthEast'); |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 | function [xplus, converged, iter, k] = butcherintegration(x, func, h, a, b, Niter) % [xplus, converged, iter, k] = butcherintegration(x, f, h, a, b, [Niter=100]) % % Performs numerical integration of f starting from x for h time units. % % f must be a function handle that returns dx/dt. If a defines an implicit % method, then f must also return the Jacobian of dx/dt. In either case, f must % be time-invariant. % % a and b should be the coefficients from the Butcher tableau for the % integration scheme. If the a matrix defines an explicit scheme (i.e., zero on % and above the diagonal, then the explicit formulas will be used. Otherwise, % the implicit equations are solved using Newton's method. % % Niter is the maximum number of iterations to perform. xguess is a guess for % x at the next timestep. narginchk(5, 6); if nargin() < 6 Niter = 100; end % Check arguments. x = x(:); Nx = length(x); b = h*b(:); Ns = length(b); if ~isequal(size(a), [Ns, Ns]) error('a must be a square matrix with a column for each element of b!'); end a = h*a; % Decide if tableau is explicit or not. if nnz(triu(a)) == 0 k = butcher_explicit(x, func, a); converged = true(); iter = 0; else [k, converged, iter] = butcher_implicit(x, func, a, Niter); end % Advance state. xplus = x + k*b; end%function % ***************************************************************************** function k = butcher_explicit(x, func, a) % Performs an explicit integration step. Nx = length(x); Ns = size(a, 1); k = zeros(Nx, Ns); for i = 1:Ns k(:,i) = func(x + k*a(i,:)'); end end%function % ***************************************************************************** function [k, converged, iter] = butcher_implicit(x, func, a, Niter) % Solves implicit integration using Newton's Method. Timestep is assumed 1. narginchk(4, 4); Nx = length(x); Ns = size(a, 1);; k = zeros(Nx, Ns); converged = false(); iter = 0; while ~converged && iter < Niter iter = iter + 1; dX = k*a'; if Ns == 1 [f, J] = func(x + dX); J = a*J; else f = zeros(Nx*Ns, 1); J = zeros(Nx*Ns, Nx*Ns); for i = 1:Ns I = ((i - 1)*Nx + 1):(i*Nx); [f(I), Ji] = func(x + dX(:,i)); J(I,:) = kron(a(i,:), Ji); end end f = f - k(:); J = J - eye(Nx*Ns); dk = -J\f; k = k + reshape(dk, [Nx, Ns]); converged = norm(dk) < 1e-8; end end%function |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 | function errorplot2d(x, T, location) % errorplot2d(x, T, [location]) % % Makes a plot of integration error over the given time interval T. % % x must be a struct with fields to plot. All step sizes are assumed to be % uniform. If x contains an 'exact' field, it is handled specially. narginchk(2, 3); if nargin() < 3 location = {}; else location = {'location', location}; end % Create figure and axes. figure(); ax1 = subplot(2, 2, 1); hold(ax1, 'on'); ylabel(ax1, 'x_1', 'rotation', 0); xlabel(ax1, 't'); ax2 = subplot(2, 2, 3); hold(ax2, 'on'); ylabel(ax2, 'x_2', 'rotation', 0); xlabel(ax2, 't'); axp = subplot(2, 2, [2, 4]); hold(axp, 'on'); ylabel(axp, 'x_2', 'rotation', 0); xlabel(axp, 'x_1'); % Plot exact solution. if isfield(x, 'exact') t = linspace(0, T, size(x.exact, 2)); plot(ax1, t, x.exact(1,:), '-k'); plot(ax2, t, x.exact(2,:), '-k'); plot(axp, x.exact(1,:), x.exact(2,:), '-k'); exact = {'exact'}; else exact = {}; end % Plot other solutions. fields = setdiff(fieldnames(x), {'exact'}); markers = {'o', 's', 'd', '^', 'v', '>', '<', 'p', 'h'}; markers = markers(mod(0:(length(fields) - 1), length(markers)) + 1); colors = hsv(length(fields)); for i = 1:length(fields) f = fields{i}; t = linspace(0, T, size(x.(f), 2)); style = {markers{i}, 'color', colors(i,:)}; plot(ax1, t, x.(f)(1,:), style{:}); plot(ax2, t, x.(f)(2,:), style{:}); plot(axp, x.(f)(1,:), x.(f)(2,:), style{:}); end % Add legend. legend(axp, exact{:}, fields{:}, location{:}); end%function |