Matlab: метод Ньютона для приближенного решения

Question

Я хотел бы реализовать приведенную ниже систему уравнений в Matlab для аппроксимации колебаний маятника. Хотя я понимаю начальную настройку, я не понимаю, как получается «дельта» при делении двух матриц!

Как деление двух матриц -G и J дает дельту (k)? Я попробовал первый шаг в Matlab, и в результате 5by5!

Как решить для дельты. Простой пример был бы великолепен.

Пример MATLAB:

theta= 0.85;
A = (1/h^2).*(diag((-2+sin(theta))*ones(1,n)) + diag(ones(1,n-1),1) + diag(ones(1,n-1),-1));
J= (1/h^2).*(diag((-2+h^2*(cos(theta))^2)*ones(1,n)) + diag(ones(1,n-1),1) + diag(ones(1,n-1),-1));

delta=-A\J

Мне просто нужно понять этот шаг, чтобы продолжить итерации.

Заранее спасибо!

Answer 1

В математической части G - это не матрица, а вектор. Точнее, это m функции, [G1(theta), G2(theta), ..., Gm(theta)], которые даются в первой строке вашей математики. Аргумент theta здесь также является вектором [theta1, ..., thetam]. Вы хотите найти m тэты, чтобы сделать вектор G всеми нулями. Вот что написано прямо перед уравнением для J.

Каждый из этих Gi имеет производные по theta1, ..., thetam. Они записаны в матричной форме как J в математике.

Чтобы решить эту нелинейную систему уравнений методом Ньютона, у вас есть предположение для некоторых значений для theta1, ..., thetam который называется theta^[k] в математике. Оцените G, чтобы получить вектор фактических чисел («остатки»). Это G(theta^[k]) в математике.

Если бы функция была линейной, вы могли бы выяснить, как точно исправить предположение, добавив вектор J(theta^[k]) \ -G(theta^[k]). Это часть обновления Ньютона. Поскольку истинная функция является нелинейной, коррекция является лишь приближением. После вычисления theta^[k+1] вы повторяете процесс, пока не будете удовлетворены качеством ответа.

(В общем случае метод Ньютона не обязательно сходится, но если вы находитесь в окрестности решение и / или вещи не слишком нелинейны, это довольно быстро, удваивая количество правильных цифр на каждой итерации.)

Редактировать: пример для m=5

m=5;
x=linspace(0, 3, m)';
theta=1/5*cos(x)+1/2*sin(x);
alpha=0;
beta=-0.2;
h=0.1;
steps=3;
format short e;
for k=1:steps+1
  theta_minus1=[alpha; theta(1:end-1)];
  theta_plus1=[theta(2:end); beta];
  G=1/h^2*(theta_minus1-2*theta+theta_plus1) + sin(theta);
  fprintf('G=');
  disp(G');
  if k == steps+1
    break
  end
  J=1/h^2*(diag(-2+h^2*cos(theta)) + diag(ones(1, m-1), 1) + diag(ones(1, m-1), -1));
  delta=J\-G;
  theta=theta+delta;
end