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 | % Illustration of GL4 collocation approximation % Butcher tableau for GL4. tab = struct('a', [1/4, 1/4 - sqrt(3)/6; 1/4 + sqrt(3)/6, 1/4], ... 'b', [1/2, 1/2], ... 'c', [1/2 - sqrt(3)/6, 1/2 + sqrt(3)/6]); % Run one step of GL4. h = 2*pi()/5; x0 = [1; 0]; [x1, ~, ~, k] = butcherintegration(x0, @stiffode, h, tab.a, tab.b); % Get exact solution. t = linspace(0, h, 1001); x = [cos(t); -sin(t)]; xdot = [-sin(t); -cos(t)]; % Get collocation polynomial approximation. k1 = k(:,1); k2 = k(:,2); c1 = tab.c(1); c2 = tab.c(2); rho1 = h*(k1*c2 - k2*c1)/(c1 - c2); rho2 = 0.5*h*(k1 - k2)/(c1 - c2); xc = NaN(2, length(t)); xdotc = NaN(2, length(t)); xnode = NaN(2, length(tab.c)); for i = 1:2 xc(i,:) = x0(i) - rho1(i)*(t/h) + rho2(i)*(t/h).^2; xdotc(i,:) = -rho1(i)/h + 2*rho2(i)/h*(t/h); xnode(i,:) = x0(i) - rho1(i)*tab.c + rho2(i)*tab.c.^2; end % Get tangent lines for derivative. dt = 0.03; % Width for tangents. tt = [tab.c(1) + [-dt, dt], NaN(), tab.c(2) + [-dt, dt]]; xt = [xnode(:,1) - h*k1*dt, xnode(:,1) + h*k1*dt, NaN(2, 1), ... xnode(:,2) - h*k2*dt, xnode(:,2) + h*k2*dt]; % Make a plot. figure(); t = t/h; for i = 1:2 subplot(2, 2, i); hold('on'); plot(t, x(i,:), '-k', t, xc(i,:), '--k'); plot(tt, xt(i,:), '-k', 'linewidth', 5); ylabel(sprintf('x_%d', i), 'rotation', 0); xlabel('\tau/h'); subplot(2, 2, 2 + i); hold('on'); plot(t, xdot(i,:), '-k', t, xdotc(i,:), '--k'); plot(tab.c, [k1(i), k2(i)], 'ok', 'markerfacecolor', 'k'); ylabel(sprintf('dx_%d/dt', i), 'rotation', 0); xlabel('\tau/h'); end % Save data. data = struct('t', t, 'x', x, 'xdot', xdot, 'xc', xc, 'xdotc', xdotc, ... 'tt', tt, 'xt', xt, 'h', h, 'k', k, 'c', tab.c, 'xn', xnode); save('-v7', 'gl4illustration.mat', '-struct', 'data'); |
1 2 3 4 5 6 7 8 9 10 | function [xdot, jac] = stiffode(x) % [xdot, jac] = stiffode(x) % % Example stiff ODE. First output is dx/dt; second is Jacobian of dx/dt. A = [0 1; -1 0]; lambda = 500; rho = (x'*x - 1); xdot = A*x - lambda*rho*x; jac = A - eye(2)*lambda*rho - 2*lambda*(x*x'); 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 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 |