Метод Гаусса – Лежандра - Gauss–Legendre method

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

Все методы Гаусса – Лежандра A-устойчивы .

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

1/2	1/2
	1

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

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

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

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

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

Интуиция

Методы Гаусса-Лежандра используются для решения любого обыкновенного дифференциального уравнения вида , где - некоторый вектор в . Интуитивная идея состоит в том, чтобы взять точное тождество и аппроксимировать интеграл в правой части квадратурой Гаусса . Например, использование двухточечной квадратуры дает приближение ${\ Displaystyle {\ точка {х}} = е (х)}$${\ Displaystyle х (т)}$${\ Displaystyle \ mathbb {R} ^ {п}}$${\ displaystyle x (h) -x (0) = \ int _ {0} ^ {h} f (x (t)) \, dt}$

${\ displaystyle x (h) -x (0) = {\ frac {1} {2}} f \ left (x \ left (\ left ({\ frac {1} {2}} - {\ frac {\ sqrt {3}} {6}} \ right) h \ right) \ right) + {\ frac {1} {2}} f \ left (x \ left (\ left ({\ frac {1} {2} } + {\ frac {\ sqrt {3}} {6}} \ right) h \ right) \ right) + O (h ^ {4})}$.

К сожалению, необходимы дальнейшие приближения, так как правую часть по-прежнему невозможно точно оценить. Мы ищем приближение Рунге – Кутты для оценки в эти моменты времени. Позвольте и быть сокращениями для этих двух квадратурных точек. ${\ displaystyle x}$${\ displaystyle x_ {1}}$${\ displaystyle x_ {2}}$

${\ displaystyle x_ {1} = x \ left (\ left ({\ frac {1} {2}} - {\ frac {\ sqrt {3}} {6}} \ right) h \ right) \ приблизительно x (0) + a_ {11} k_ {1} + a_ {12} k_ {2}.}$

${\ displaystyle x_ {2} = x \ left (\ left ({\ frac {1} {2}} + {\ frac {\ sqrt {3}} {6}} \ right) h \ right) \ приблизительно x (0) + a_ {21} k_ {1} + a_ {22} k_ {2}.}$

Вот . Это неявное ограничение, которое должно быть решено с помощью алгоритма поиска корня, такого как метод Ньютона . Значения параметров Рунге-Кутты можно определить из разложения в ряд Тейлора в . Полученная схема является методом четвертого порядка. ${\ Displaystyle к_ {я} = е (х_ {я})}$${\ displaystyle a_ {ij}}$${\ displaystyle h}$

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

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

%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

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

Интегрированная орбита вблизи аттрактора Лоренца.

Примечания

использованная литература

• Изерлес, Ари (1996), Первый курс численного анализа дифференциальных уравнений , Cambridge University Press , ISBN 978-0-521-55655-2.