Метод Гаусса – Лежандра

В численном анализе и научных вычислениях методы Гаусса -Лежандра представляют собой семейство численных методов для обыкновенных дифференциальных уравнений . Методы Гаусса–Лежандра являются неявными методами Рунге–Кутты . Точнее, это методы коллокации, основанные на точках квадратуры Гаусса – Лежандра . Метод Гаусса–Лежандра, основанный на s точках, имеет порядок 2 s . ^{[ 1 ]}

Все методы Гаусса–Лежандра A-стабильны . ^{[ 2 ]}

Метод Гаусса-Лежандра второго порядка представляет собой неявное правило средней точки . Его таблица Мясника :

1/2	1/2
	1

Метод Гаусса – Лежандра четвертого порядка имеет таблицу Мясника:

${\tfrac {1}{2}}-{\tfrac {1}{6}}{\sqrt {3}}$	${\tfrac {1}{4}}$	${\tfrac {1}{4}}-{\tfrac {1}{6}}{\sqrt {3}}$
${\tfrac {1}{2}}+{\tfrac {1}{6}}{\sqrt {3}}$	${\tfrac {1}{4}}+{\tfrac {1}{6}}{\sqrt {3}}$	${\tfrac {1}{4}}$
	${\tfrac {1}{2}}$	${\tfrac {1}{2}}$

Метод Гаусса – Лежандра шестого порядка имеет таблицу Мясника:

${\tfrac {1}{2}}-{\tfrac {1}{10}}{\sqrt {15}}$	${\tfrac {5}{36}}$	${\tfrac {2}{9}}-{\tfrac {1}{15}}{\sqrt {15}}$	${\tfrac {5}{36}}-{\tfrac {1}{30}}{\sqrt {15}}$
${\tfrac {1}{2}}$	${\tfrac {5}{36}}+{\tfrac {1}{24}}{\sqrt {15}}$	${\tfrac {2}{9}}$	${\tfrac {5}{36}}-{\tfrac {1}{24}}{\sqrt {15}}$
${\tfrac {1}{2}}+{\tfrac {1}{10}}{\sqrt {15}}$	${\tfrac {5}{36}}+{\tfrac {1}{30}}{\sqrt {15}}$	${\tfrac {2}{9}}+{\tfrac {1}{15}}{\sqrt {15}}$	${\tfrac {5}{36}}$
	${\tfrac {5}{18}}$	${\tfrac {4}{9}}$	${\tfrac {5}{18}}$

Вычислительные затраты на методы Гаусса – Лежандра более высокого порядка обычно чрезмерны, и поэтому они используются редко. ^{[ 3 ]}

Интуиция

Методы Гаусса-Лежандра, Рунге-Кутты (GLRK) решают обыкновенное дифференциальное уравнение. ${\dot {x}}=f(x)$ с $x(0)=x_{0}$ . Отличительной особенностью GLRK является оценка ${\textstyle x(h)-x_{0}=\int _{0}^{h}dt\,f(x(t))}$ с гауссовой квадратурой .

$x(h)=x(0)+{\frac {h}{2}}\sum _{i=1}^{\ell }w_{i}k_{i}+O(h^{2\ell }),$

где $k_{i}=f(x(hc_{i}))$ являются выборочными скоростями, $w_{i}$ квадратурные веса, ${\textstyle c_{i}={\frac {1}{2}}h(1+r_{i})}$ представляют собой абсциссы, а $r_{i}$ это корни $P_{\ell }(r_{i})=0$ полинома Лежандра степени $\ell$ . Необходимо дальнейшее приближение, поскольку $k_{i}$ пока оценить невозможно. Чтобы сохранить ошибку усечения порядка $O(h^{2\ell })$ , нам нужно только $k_{i}$ заказать $O(h^{2\ell -1})$ . Неявное определение Рунге-Кутты ${\textstyle k_{i}=f{\left(x_{0}+h\sum _{j}a_{ij}k_{j}\right)}}$ вызывается для достижения этой цели. Это неявное ограничение, которое необходимо решить с помощью алгоритма поиска корня, такого как метод Ньютона . Значения параметров Рунге-Кутты $a_{ij}$ можно определить из разложения в ряд Тейлора по $h$ .

Практический пример

Методы Гаусса-Лежандра неявны, поэтому в общем случае их нельзя точно применить. Вместо этого делается обоснованное предположение о $k_{i}$ , а затем использует метод Ньютона для сходимости сколь угодно близко к истинному решению. Ниже приведена функция Matlab, реализующая метод Гаусса-Лежандра четвертого порядка.

% starting point
x = [ 10.5440; 4.1124; 35.8233];

dt = 0.01;
N = 10000;
x_series = [x];
for i = 1:N
  x = gauss_step(x, @lorenz_dynamics, dt, 1e-7, 1, 100);
  x_series = [x_series x];
end

plot3( x_series(1,:), x_series(2,:), x_series(3,:) );
set(gca,'xtick',[],'ytick',[],'ztick',[]);
title('Lorenz Attractor');
return;

function [td, j] = lorenz_dynamics(state)
  % return a time derivative and a Jacobian of that time derivative
  x = state(1);
  y = state(2);
  z = state(3);

  sigma = 10;
  beta  = 8/3;
  rho   = 28;

  td = [sigma*(y-x); x*(rho-z)-y; x*y-beta*z];

  j = [-sigma, sigma, 0;
        rho-z, -1, -x;
        y, x, -beta];
end

function x_next = gauss_step( x, dynamics, dt, threshold, damping, max_iterations )
  [d,~] = size(x);
  sq3 = sqrt(3);
  if damping > 1 || damping <= 0
    error('damping should be between 0 and 1.')
  end

  % Use explicit Euler steps as initial guesses
  [k,~] = dynamics(x);
  x1_guess = x + (1/2-sq3/6)*dt*k;
  x2_guess = x + (1/2+sq3/6)*dt*k;
  [k1,~] = dynamics(x1_guess);
  [k2,~] = dynamics(x2_guess);

  a11 = 1/4;
  a12 = 1/4 - sq3/6;
  a21 = 1/4 + sq3/6;
  a22 = 1/4;

  error = @(k1, k2) [k1 - dynamics(x+(a11*k1+a12*k2)*dt); k2 - dynamics(x+(a21*k1+a22*k2)*dt)];
  er = error(k1, k2);
  iteration=1;
  while (norm(er) > threshold && iteration < max_iterations)
    fprintf('Newton iteration %d: error is %f.\n', iteration, norm(er) );
    iteration = iteration + 1;
   
    [~, j1] = dynamics(x+(a11*k1+a12*k2)*dt);
    [~, j2] = dynamics(x+(a21*k1+a22*k2)*dt);
   
    j = [eye(d) - dt*a11*j1, -dt*a12*j1;
         -dt*a21*j2, eye(d) - dt*a22*j2];
    
    k_next = [k1;k2] - damping * linsolve(j, er);
    k1 = k_next(1:d);
    k2 = k_next(d+(1:d));
   
    er = error(k1, k2);
  end
  if norm(er) > threshold
    error('Newton did not converge by %d iterations.', max_iterations);
  end
  x_next = x + dt / 2 * (k1 + k2);
end

Этот алгоритм на удивление дешев. Ошибка в $k_{i}$ может упасть ниже $10^{-12}$ всего за 2 шага Ньютона. Единственная дополнительная работа по сравнению с явными методами Рунге-Кутты — это вычисление якобиана.

Симметричные во времени варианты

За счет добавления дополнительного неявного соотношения эти методы можно адаптировать для обеспечения симметрии обращения времени. В этих методах усредненное положение $(x_{f}+x_{i})/2$ используется в вычислениях $k_{i}$ вместо просто начальной позиции $x_{i}$ стандартными методами Рунге-Кутты. Метод порядка 2 — это всего лишь неявный метод средней точки.