Soustava nelineárních rovnic F(X) = 0

Newtonova metoda

Příklad: hledáme řešení soustavy rovnic
x2 = x13 + 1
x12 + x22 = 4

Soustavu napíšeme ve tvaru F(X) = 0, kde X = [ x1 , x2 ]T, F = [ f1(X) , f2(X) ]T,
f1(X) = x12 + x22 - 4 ,
f2(X) = x2 - x13 - 1 .

V Matlabu zadáme F jako

>> syms x1 x2          % definujeme promenne x1, x2
>> X = [ x1 ; x2]      % vektor promennych x1, x2
>> f1 = x1^2 + x2^2 -4 % funkce f1 promennych x1, x2
>> f2 = x2-x1^3-1      % funkce f2 promennych x1, x2
>> F = [ f1 ; f2 ]     % vektorova funkce F promennych x1, x2

Pro Newtonovu metodu nejdřív potřebujeme spočítat Jacobiovu matici funkce F

>> J=jacobian(F, X)        % takhle
>> J=jacobian(F, [x1, x2]) % nebo takhle (pak nemusime vytvaret vektor X)

a zvolit počáteční aproximaci X0:

>> X0 = [1; 2]

Pak opakujeme následující 4 kroky:

>> J0 = subs( J, {x1, x2}, {X0(1), X0(2)}) % do J dosadime za X hodnotu X0
>> F0 = subs(F, {x1, x2}, {X0(1), X0(2)})  % do F dosadime za X hodnotu X0
>> D0 = J0 \ (- F0)                        % resime soustavu J0 * D0 = - F0
>> X0 = X0 + D0                            % opravime aproximaci X0


Prostá iterační metoda (PIM)

Řešíme stejnou úlohu jako výše, případně pozměněnou úlohu, kde první rovnici nahradíme rovnicí
x2 = sin(x1) + 1

V následujícím skriptu si zkuste měnit hodnoty proměnných pit (počet iterací), X (výchozí iterace) a alfa (pomocný parametr, kterým ovlivňujeme konvergenci).

syms x1 x2 ;         % definujeme promenne x1, x2
X = [ x1 ; x2];      % vektor promennych x1, x2
f1 = x1^2 + x2^2 -4; % funkce f1 promennych x1, x2
f2 = x2-x1^3-1;      % funkce f2 promennych x1, x2
% f2 = x2-sin(x1)-1 ; % jina uloha; zde je lepsi alfa = 0.25, viz nize
F = [ f1 ; f2 ]     % vektorova funkce F promennych x1, x2

pit = 30;  % max. pocet iteraci
X = [1; 2]; % vychozi aproximace
alfa = 0.1;    % koeficient pro lepsi konvergenci - hadame
G = [x1;x2] - alfa*F; % budeme resit rovnici x = G(x), tj. F(x) = 0
for k=1:maxit
 X = subs(G, {x1, x2}, {X(1), X(2)}); % x_{k} = G(x_{k-1})
 % zobrazeni poslednich 3 iteraci
 if (k>(pit-3))
   if(k==(pit-2))
     disp('posledni 3 aproximace X:')
   end
   disp(X)
 end
end

% zkouska
disp('hodnota F(X):') 
FX = subs(F, {x1, x2}, {X(1), X(2)});
disp(FX);