Correzione codice C++ con newton raphson (matlab)

[-Root-]

ragazzi, c'è qualche buon'anima che mi può correggere dove sbaglio con questo codice per risolvere un sistema di equazioni non lineari con newton raphson?

[function [r,F,niter] = NR1(x0, tol, maxiter)
% Metodo di Newton-Raphson per sistemi di equazioni non lineari
% r = zero del sistema di equazioni non lineari f(x) = 0
% x0 = approssimazione iniziale dello zero r
% Il calcolo viene interrotto sia se la norma di f
28
% all'approsimazione corrente risulta essere minore di tol
% sia se l'errore relativo tra due soluzioni successive
% risulta essere minore di tol, o se viene raggiunto
% il numero massimo di iterazioni maxiter.
 % reciproco del numero di condizionamento
% Se J è mal condizionata rcond è vicino ad eps
tol = 0.1;
xnew = x(0);
if (norm(fun0(xnew))<tol)&(norm(xnew(k)-xnew,'inf')/norm(xnew1,'inf')<tol);
   xnew(k) = xnew;
else
    for k=1:maxiter
        xnew(k) = xnew(k-1) - Ji(xnew(k-1))*fun0(xnew(k-1)) % risolvo il sistema lineare
        if (norm(fun0(xnew(k))<tol)&(norm(xnew(k)-xnew(k-1)),'inf')/norm(xnew(k),'inf')<tol);
     break
        end
    end
niter = k; % numero di iterazioni effettivamente eseguite
r = xnew(k); % zero del sistema di equazioni non lineari
F = fun0(xnew1);
end]

[Daniels118]

Dove hai valorizzato x(0)?
Se la prima condizione è soddisfatta quanto vale k?
A che serve salvare in un vettore i risultati di tutte le iterazioni?
Quanto vale xnew(k-1) alla prima iterazione?
Dove valorizzi xnew1?

[-Root-]

Originariamente inviato da Daniels118

Dove hai valorizzato x(0)?
Se la prima condizione è soddisfatta quanto vale k?
A che serve salvare in un vettore i risultati di tutte le iterazioni?
Quanto vale xnew(k-1) alla prima iterazione?
Dove valorizzi xnew1?

allora, in sostanza devo valutare se il sistema in un punto x0 = [..;..] scelto da me, sia inferiore a un Epsilon sempre scelto da me, se si allora mi restituisce il risultato, se no passo alla prossima variabile trovata con x=x0 - F/J dove J è lo jacobiano (matrice 2 x 2), e F è il vettore colonna delle due funzioni che compongono il sistema. E così via fino a che x(k) = x(k-1) - F/J e se F(x(K)) è minore di epsilon, e la norma di F/J è minore di un delta sempre scelto da me allora mi da il valore di x(k) e della funzione F(x(k)). Quel codice di prima non lo sto usando più, ne ho fatto uno nuovo, che però quando lo faccio partire (inserisco nmax (n max iterazioni), tol (sarebbe la epsilon, che per comodità l'ho resa uguale a delta anche, e x0, e poi scrivo la funzione) mi dice Matrix dimensions must agree. (relativo a F/J) . questo è il codice

[function [x,res,niter,Err] = newtonsys(Ffun,Jfun,x0,tol,...
                                nmax, varargin)
%NEWTONSYS calcola una radice di un sistema non lineare
%   [ZERO,RES,NITER,ERR]=NEWTONSYS(FFUN,JFUN,X0,TOL,NMAX)
%   calcola il vettore ZERO, radice di un sistema non
%   lineare definito nella function FFUN con matrice
%   jacobiana definita nella function JFUN a partire
%   dal vettore X0. RES contiene il valore del residuo
%   in ZERO e NITER il numero di iterazioni necessarie
%   per calcolare ZERO. FFUN e JFUN sono function
%   definite tramite M-file
%   Nel vettore ERR sono salvate le norme degli incrementi
%   alle varie iterazioni

niter = 0;
Func = tol + 1;
err = tol + 1;
x = x0;
Err=[];
while Func >= tol & err >= tol & niter < nmax
    J = Jfun(x,varargin{:});
    F = Ffun(x,varargin{:});
    z = - J/F;
    x = x + z;
    err = norm(z);
    Func = norm(F);
    Err=[Err;err];
    niter = niter + 1;
end
res = norm(Ffun(x,varargin{:}));
if (niter==nmax & err> tol)
    fprintf(['Il metodo non converge nel massimo ',...
       'numero di iterazioni. L''ultima iterata\n',...
       'calcolata ha residuo relativo pari a %e\n'],res);
else
    fprintf(['Il metodo converge in %i iterazioni',...
            ' con un residuo pari a %e\n'],niter,res);
end
return]

ps queste sono le function del sistema e della matrice jacobiana

[function F = Ffun(x)
  F(1,1) = 0.3*cos(0) + 0.5*cos(x(1)) - 0.7*cos(x(2)) - 0.4;
  F(2,1) = 0.3*sin(0) + 0.5*sin(x(1)) + 0.7*sin(x(2));
return]  [function J = Jfun(x)
J(1,1)= -0.5*sin(x(1))
J(1,2)= -0.7*sin(x(2))
J(2,1)= 0.5*cos(x(1))
J(2,2)= 0.7*sin(x(2))
return]

[-Root-]

Originariamente inviato da Antonio23

ma perché il codice è sotto spoiler ?! devi usare i tag [ code ]

non avevo visto che c'era il comando code chiedo venia e riparo subito

[Daniels118]

Per favore quando posti il codice usa i tag code, così è illeggibile... per questa volta lo riscrivo io:

[function [x,res,niter,Err] = newtonsys(Ffun,Jfun,x0,tol,nmax, varargin)
%NEWTONSYS calcola una radice di un sistema non lineare
% [ZERO,RES,NITER,ERR]=NEWTONSYS(FFUN,JFUN,X0,TOL,NMAX)
% calcola il vettore ZERO, radice di un sistema non
% lineare definito nella function FFUN con matrice
% jacobiana definita nella function JFUN a partire
% dal vettore X0. RES contiene il valore del residuo
% in ZERO e NITER il numero di iterazioni necessarie
% per calcolare ZERO. FFUN e JFUN sono function
% definite tramite M-file
% Nel vettore ERR sono salvate le norme degli incrementi
% alle varie iterazioni
niter = 0;
Func = tol + 1;
err = tol + 1;
x = x0;
Err=[];
while Func >= tol & err >= tol & niter < nmax
  J = Jfun(x,varargin{:});
  F = Ffun(x,varargin{:});
  z = - J/F;
  x = x + z;
  err = norm(z);
  Func = norm(F);
  Err=[Err;err];
  niter = niter + 1;
end
res = norm(Ffun(x,varargin{:}));
if (niter==nmax & err> tol)
  fprintf(['Il metodo non converge nel massimo ',...
  'numero di iterazioni. L''ultima iterata\n',...
  'calcolata ha residuo relativo pari a %e\n'],res);
else
  fprintf(['Il metodo converge in %i iterazioni',...
  ' con un residuo pari a %e\n'],niter,res);
end
return]

[function F = Ffun(x)
  F(1,1) = 0.3*cos(0) + 0.5*cos(x(1)) - 0.7*cos(x(2)) - 0.4;
  F(2,1) = 0.3*sin(0) + 0.5*sin(x(1)) + 0.7*sin(x(2));
return]

[function J = Jfun(x)
  J(1,1)= -0.5*sin(x(1))
  J(1,2)= -0.7*sin(x(2))
  J(2,1)= 0.5*cos(x(1))
  J(2,2)= 0.7*sin(x(2))
return]

Prima hai scritto di dover dividere F per J, ma nel codice c'è scritto J/F... qual è la strada corretta?

[-Root-]

Originariamente inviato da Daniels118

Per favore quando posti il codice usa i tag code, così è illeggibile... per questa volta lo riscrivo io:

[function [x,res,niter,Err] = newtonsys(Ffun,Jfun,x0,tol,nmax, varargin)
%NEWTONSYS calcola una radice di un sistema non lineare
% [ZERO,RES,NITER,ERR]=NEWTONSYS(FFUN,JFUN,X0,TOL,NMAX)
% calcola il vettore ZERO, radice di un sistema non
% lineare definito nella function FFUN con matrice
% jacobiana definita nella function JFUN a partire
% dal vettore X0. RES contiene il valore del residuo
% in ZERO e NITER il numero di iterazioni necessarie
% per calcolare ZERO. FFUN e JFUN sono function
% definite tramite M-file
% Nel vettore ERR sono salvate le norme degli incrementi
% alle varie iterazioni
niter = 0;
Func = tol + 1;
err = tol + 1;
x = x0;
Err=[];
while Func >= tol & err >= tol & niter < nmax
  J = Jfun(x,varargin{:});
  F = Ffun(x,varargin{:});
  z = - J/F;
  x = x + z;
  err = norm(z);
  Func = norm(F);
  Err=[Err;err];
  niter = niter + 1;
end
res = norm(Ffun(x,varargin{:}));
if (niter==nmax & err> tol)
  fprintf(['Il metodo non converge nel massimo ',...
  'numero di iterazioni. L''ultima iterata\n',...
  'calcolata ha residuo relativo pari a %e\n'],res);
else
  fprintf(['Il metodo converge in %i iterazioni',...
  ' con un residuo pari a %e\n'],niter,res);
end
return]

[function F = Ffun(x)
  F(1,1) = 0.3*cos(0) + 0.5*cos(x(1)) - 0.7*cos(x(2)) - 0.4;
  F(2,1) = 0.3*sin(0) + 0.5*sin(x(1)) + 0.7*sin(x(2));
return]

[function J = Jfun(x)
  J(1,1)= -0.5*sin(x(1))
  J(1,2)= -0.7*sin(x(2))
  J(2,1)= 0.5*cos(x(1))
  J(2,2)= 0.7*sin(x(2))
return]

Prima hai scritto di dover dividere F per J, ma nel codice c'è scritto J/F... qual è la strada corretta?

ho controllato meglio il codice da cui ho preso spunto, li era scritto -J\F, anzi che /, in teoria dovrebbe essere la matrice colonna F moltiplicata per l0inversa dello jacobiano. Ora il programma funziona, non mi da messaggi di errore, ma mi da sempre il massimo numero di iterazioni, e la funzione non è mai minore dell'epsilon che ho scelto, inoltre vengono soluzioni troppo grandi, tenendo presente che x(2) e x(1) sono due angoli, (per info, è l'analisi cinematica di un quadrilatero articolato, un meccanismo) . Ho provato a cambiare il punto iniziale ma niente, ora il codice è così:

function [x,res,niter,Err] = newtonsys(Ffun,Jfun,x0,tol,...
                                nmax, varargin)
%NEWTONSYS calcola una radice di un sistema non lineare
%   [ZERO,RES,NITER,ERR]=NEWTONSYS(FFUN,JFUN,X0,TOL,NMAX)
%   calcola il vettore ZERO, radice di un sistema non
%   lineare definito nella function FFUN con matrice
%   jacobiana definita nella function JFUN a partire
%   dal vettore X0. RES contiene il valore del residuo
%   in ZERO e NITER il numero di iterazioni necessarie
%   per calcolare ZERO. FFUN e JFUN sono function
%   definite tramite M-file
%   Nel vettore ERR sono salvate le norme degli incrementi
%   alle varie iterazioni

niter = 0;
Func = tol;
Jac = tol;
x = x0;
Err=[];
while Func >= tol  & niter < nmax & Jac >= tol
    J = Jfun(x,varargin{:});
    F = Ffun(x,varargin{:});
    z = - J\F;
    x = x + z;
    Jac = norm(J);
    Func = norm(F);
    niter = niter + 1;
end
res = norm(Ffun(x,varargin{:}));
if (Func > tol)
    fprintf(['Il metodo non converge nel massimo ',...
       'numero di iterazioni. L''ultima iterata\n',...
       'calcolata ha residuo relativo pari a %e\n'],res);
else
    fprintf(['Il metodo converge in %i iterazioni',...
            ' con un residuo pari a %e\n'],niter,res);
end
return