Как я итерирую / увеличиваю, чтобы найти корни для системы уравнений с двумя переменными, используя метод Ньютона - PullRequest
Новая
       17

Как я итерирую / увеличиваю, чтобы найти корни для системы уравнений с двумя переменными, используя метод Ньютона

0 голосов
/ 10 апреля 2019

Как я уже говорил в заголовке, мне интересно найти корень для системы из двух уравнений и двух переменных, а именно:

f (x, y) = 0
g (x, y)= 0

для функций:

f (x, y) = x (1 + y ^ 2) -1
g (x, y) = y (1 + x^ 2) -2

Как правильно выполнить итерацию, так как я подозреваю, что мой код в этом отношении разбивается на части?

Я подтвердил, вычислив это вручную в командной строке и командойСудя по всему, я получаю довольно хорошие приближения для корня после 4 итераций.Когда допуск достаточно близок к 0.

Кроме этого, я только что изменил длину вектора и различные значения.Ничего не стоит упоминать.

function [x,y] = NewtonMeth(f, g, a, b)

n=100; %arbitrary amount of iterations
x=zeros(1,n+1); %initializing my vector for x values 
y=zeros(1,n+1);  %initializing my vector for y values
y(1)=b; %setting my first element as the intial guess
x(1)=a;


tol = 0.00001; %tolerance, denoting when to stop my for loop.

for i=1:n
%utilizing a previous function for partial derivative:

   f1=PD(f, x(i), y(i), 1); %partial derivatives for my two functions f & g 
   f2=PD(f, x(i), y(i), 2);
   g1=PD(g, x(i), y(i), 1);
   g2=PD(g, x(i), y(i), 2);

   j1= [f(x(i),y(i)) f2; g(x(i),y(i)) g2]; %separating my jacobian into amatrix containing the relevant elements
   j2 = [f1 f2; g1 g2];%second jacobian matrix
   j = det(j1)./det(j2); %here I evaluate the jacobian, with accordance to the formula for the newton method for system of equations 

   %with two variables. 

   if abs(f(x(i),y(i)))>=tol %Said tolerance which should limit the amount of iterations

      x(i+1)= x(i)-j; %incrementing the x values inorder to approximate the root 

      y(i+1)=y(i)-j; %incrementing the y values inorder to approximate the root

   end
end
end

Вот моя функция для получения частной производной: 

function derivative = PD(f, a, b, i)
h = 0.000001;
fn=zeros(1,2);
if i == 1
    fn(i) = (f(a+h,b)-f(a,b))/h;
elseif i==2
    fn(i) = (f(a,b+h)-f(a,b))/h;
end
derivative = fn(i);
end

Выбор начального предположения для a=0.2 и b=1.8 (как сделано в моей книге Calc 3) Iдолжен получить: 

|x(1) = 0.200000|y(1)=1.800000|f(x(1),y(1))=-0.152000|g(x(1),y(1))=-0.128000|
|x(2) = 0.216941|y(2)=1.911349|f(x(2),y(2))= 0.009481|g(x(2),y(2))= 0.001303|
|x(3) = 0.214827|y(3)=1.911779|f(x(3),y(3))=-0.000003|g(x(3),y(3))= 0.000008|
|x(4) = 0.214829|y(4)=1.911769|f(x(4),y(4))= 0.000000|g(x(4),y(4))= 0.000000|
...