[MATLAB] Sistemi non lineari con Matlab

[aeneas2019]

Dunque, devo risolvere un sistema non lineare di equazioni con il metodo di Newton-Raphson al variare di un parametro t (tempo); ho pensato di scrivere un metodo che disponga i vettori delle soluzioni come colonne di una matrice, per poi tracciarne i grafici con il comando plot.
Il listato è questo:

function T = newra(x0,eps,max_iter) 
T=zeros(61,8); 
i=1;
 t=0;
 f=fun(x0,t); 
norma_residui=norm(f);
 while(t<=6)&&(i<=61) 
k=0;
 while(norma_residui>eps)
 f=fun(x0,t);  
J=jac(x0,t); % Calcola la norma del vettore dei residui norma_residui=norm(f); 
 x=x0-J\f;
k=k+1; 
 if(k>max_iter) 
disp('no');
 disp(k); 
break, 
end 
 x0=x;
 end 
t=t+0.1;
 for j=1:8
 T(i,j)=x(j); 
end 
i=i+1; 
x0=x;
 end 
end

dove fun e jac sono metodi che sostituiscono rispettivamente al sistema di equazioni e alla matrice Jacobiana gli elementi del vettore x0 e il valore del tempo x.
Il risultato che ottengo è una matrice 61X8 (come deve essere) in cui però compare in tutte le colonne la stessa soluzione, quella con t=0. Come mai?

[rsdpzed]

non sono un utilizzatore di matlab ma cosi ad occhio interpretando il comportamento che descrivi mi viene da pensare che l'inizializzazione di t (che sembra un float visto che lo incrementi poi di 0.1) dovrebbe essere:
t = 0.0;

[aeneas2019]

ho provato ma purtroppo non cambia niente

[rsdpzed]

Originariamente inviato da aeneas2019
ho provato ma purtroppo non cambia niente

come non detto, il mio era giusto un sospetto. dovresti sentire qualcuno che conosce matlab per davvero

[drudox]

x=x0-J/f; non \

se poi posti le function e in programma main lo posso provare

inoltre a me sembra che

f=fun(x0,t);

vada dentro il while

[aeneas2019]

function f=fun(x,t)

% x: Vettore delle variabili
% f: Vettore dei residui
f(1)= 1.993*cos(x(1))+0.437*cos(x(2))+2.019*cos(x(3))+0.227*cos(0.8203);
f(2)= 1.993*sin(x(1))+0.437*sin(x(2))+2.019*sin(x(3))+0.227*sin(0.8203);
f(3)= (1.797+0.1*sin(2*pi*t/6))*cos(x(4))+0.584*cos(x(2)+0.3665)+2.019*cos(x(3))+0.908*cos(4.1713);
f(4)= (1.797+0.1*sin(2*pi*t/6))*sin(x(4))+0.584*sin(x(2)+0.3665)+2.019*sin(x(3))+0.908*sin(4.1713);
f(5)= 2.809*cos(x(1)+0.0873)+0.769*cos(x(5))+(1.54+0.1*sin(2*pi*t/6))*cos(x(6))+0.732*cos(x(1)-2.6878);
f(6)= 2.809*sin(x(1)+0.0873)+0.769*sin(x(5))+(1.54+0.1*sin(2*pi*t/6))*sin(x(6))+0.732*sin(x(1)-2.6878);
f(7)= 2.809*cos(x(1)+0.0873)-1.993*cos(x(1))-0.424*cos(x(2)+2.2340)+1.888*cos(x(5)-0.1571)+0.529*cos(x(7))+(1.907+0.1*sin(2*pi*t/6))*cos(x(8));
f(8)= 2.809*sin(x(1)+0.0873)-1.993*sin(x(1))-0.424*sin(x(2)+2.2340)+1.888*sin(x(5)-0.1571)+0.529*sin(x(7))+(1.907+0.1*sin(2*pi*t/6))*sin(x(8));

% Forma il vettore colonna

f=[f(1);f(2);f(3);f(4);f(5);f(6);f(7);f(8)];
return
end

function J=jac(x,t)

J(1,1)=-1.993*sin(x(1));
J(1,2)=-0.437*sin(x(2));
J(1,3)=-2.019*sin(x(3));
J(1,4)=0;
J(1,5)=0;
J(1,6)=0;
J(1,7)=0;
J(1,8)=0;

J(2,1)=1.993*cos(x(1));
J(2,2)=0.437*cos(x(2));
J(2,3)=2.019*cos(x(3));
J(2,4)=0;
J(2,5)=0;
J(2,6)=0;
J(2,7)=0;
J(2,8)=0;

J(3,1)=0;
J(3,2)=-0.584*sin(x(2)+0.3665);
J(3,3)=-2.019*sin(x(3));
J(3,4)=-(1.797+0.1*sin(2*pi*t/6))*sin(x(4));
J(3,5)=0;
J(3,6)=0;
J(3,7)=0;
J(3,8)=0;

J(4,1)=0;
J(4,2)=0.584*cos(x(2)+0.3665);
J(4,3)=2.019*cos(x(3));
J(4,4)=(1.797+0.1*sin(2*pi*t/6))*cos(x(4));
J(4,5)=0;
J(4,6)=0;
J(4,7)=0;
J(4,8)=0;

J(5,1)=-2.809*sin(x(1)+0.0873)-0.732*sin(x(1)-2.6878);
J(5,2)=0;
J(5,3)=0;
J(5,4)=0;
J(5,5)=-0.769*sin(x(5));
J(5,6)=-(1.54+0.1*sin(2*pi*t/6))*sin(x(6));
J(5,7)=0;
J(5,8)=0;

J(6,1)=2.809*cos(x(1)+0.0873)+0.732*cos(x(1)-2.6878);
J(6,2)=0;
J(6,3)=0;
J(6,4)=0;
J(6,5)=0.769*cos(x(5));
J(6,6)=(1.54+0.1*sin(2*pi*t/6))*cos(x(6));
J(6,7)=0;
J(6,8)=0;

J(7,1)=-2.809*sin(x(1)+0.0873)+1.993*sin(x(1));
J(7,2)=0.424*sin(x(2)+2.2340);
J(7,3)=0;
J(7,4)=0;
J(7,5)=-1.888*sin(x(5)-0.1571);
J(7,6)=0;
J(7,7)=-0.529*sin(x(7));
J(7,8)=-(1.907+0.1*sin(2*pi*t/6))*sin(x(8));

J(8,1)=2.809*cos(x(1)+0.0873)-1.993*cos(x(1));
J(8,2)=-0.424*cos(x(2)+2.2340);
J(8,3)=0;
J(8,4)=0;
J(8,5)=1.888*cos(x(5)-0.1571);
J(8,6)=0;
J(8,7)=0.529*cos(x(7));
J(8,8)=(1.907+0.1*sin(2*pi*t/6))*cos(x(8));

return
end

per attivarlo:
newra( x0, 0.00001, 10000)
dove x0=[2.81; 1.2741; 5.62; 2.1817; 5.3407; 6.1785; 5.8992; 2.0595]

Ho usato \ e non /, perché \ calcola automaticamente l'inversa della matrice e la moltiplica per il vettore; è equivalente a inv(J)*f, no?
Comunque grazie per l'aiuto


f=fun(x0, t) l'ho messo sia fuori che dentro il while; l'ho messo anche fuori per calcolarne la norma.

[drudox]

ahh si si vero ci va \ ora vedo se riesco a trovar qualcosa che non va` !

[Who am I]

A parte l' errore che non riesco a capire quale sia, dovresti fare il caso generico cioè non considerare per assunto che il vettore abbia dimensione 8, ma vedere la sua dimensione con size e operare quindi su una qualsiasi dimensione del vettore.

[aeneas2019]

Ho messo 8 perché il metodo mi serve a risolvere un problema specifico..., modificare quella cosa può fare la differenza?