trapmeth.m

来自「用matlab写的一些数值算法」· M 代码 · 共 19 行

function  [x, y] = trapmeth(A, ab, y0, N)
% Compute approximation of the solution of the initial value 
% problem  y' = Ay, y(a) = y0  on the interval ab = [a,b].
% A is a square matrix.
% Trapezoidal method with step length  h = (b-a)/N .

% Version 11.12.2003.  INCBOX

a = ab(1);  b = ab(2);  
h = (b - a)/N;  
x = linspace(a,b,N+1);     % grid
y = repmat(y0(:), 1,N+1);  % initialize  y
I = eye(size(A));          % Unit matrix
[L,U] = lu(I - h/2*A);     % LU factorization
B = I + h/2*A;             % Constant matrix on right hand side
for  n = 1 : N
  y(:,n+1) = U \ (L \ (B*y(:,n)));
end
x = x';  y = y';          % return in standard format