Sistemi lineari con MATLAB (metodo iterativo)

[andrea_d11]

Buonasera a tutti, sono nuovo e visto che sul web non ho trovato molto ho deciso di rivolgermi a voi
Da qualche giorno sono bloccato sul seguente problema:
ho un sistema lineare di 4 equazioni in 4 incognite, e fin qui tutto semplice se non fosse che le incognite sono angoli argomenti di funzioni trigonometriche, per la risoluzione quindi ho usato il metodo iterativo delle tangenti o di newton-raphson.
Sono riuscito a scrivere un programma che mi da le soluzioni (che chiamerò B e C) delle prime due equazioni e queste sono corrette. Il problema sta nel fatto che le soluzioni vengono fuori in base ad un valore che io assegno (un altro angolo chiamato A), dunque per ogni valore di A mi vengono fuori un valore di B ed uno di C; tuttavia non riesco a dire al programma di prendere e leggere i valori di A (che varia tra 0 e 360) e per ognuno darmi il valore di B e C in modo da ottenere in output un grafico di B e C in funzione di A (che è lo scopo di tutto).
Ho strutturato il programma in quattro moduli:
1.fun dove dichiaro le variabili e le equazioni
2.jacobian dove scrivo lo jacobiano associato al sistema
3.il metodo di newton raphson
4.il main (che secondo me è quello in cui c'è il problema)
Potete aiutarmi?
Grazie mille a tutti!

modulo fun

function f=fun(x)
% 
% INPUT
% x: Vettore delle incognite
% OUTPUT
% f: Vettore dei residui

r(2)=1;
r(3)=3;
r(4)=4;
r(6)=4;
r(7)=3;
for i=1:1:360   %qui ho sicuramente sbagliato qualcosa perchè in output stampa i 360 valori ti theta
    theta(1)= i;
end

f(1)= r(2)*cosd( theta(1)) + r(3)*cosd( x(1)) + r(4)*cosd( x(2)) - r(6);
f(2)= r(2)*sind( theta(1)) + r(3)*sind( x(1)) + r(4)*sind( x(2)) + r(7);
f= f';

end

modulo jacobiano

function J = jacobian(x)
%
% INPUT
% x: Vettore delle incognite
% 
% OUTPUT
% J: Matrice Jacobiana
r(3)=3;
r(4)=4;
r(2)=1;

J(1,1) = -r(3)*sind( x(1));
J(1,2) = -r(4)*sind( x(2));
J(2,1) =  r(3)*cosd( x(1));
J(2,2) =  r(4)*cosd( x(2));

end

modulo newton-raphson

function [x,iterations,f,resid,exitflag]=newtonraphson(fun,jacobian,x0,tol,itmax)
%
% INPUT
% fun : Function in cui vengono valutati gli elementi del vettore dei
% residui(equazioni di chiusura)
% jacobian : Function in cui vengono valutati gli elementi della matrice jacobiana
% x0: Vettore soluzione di primo tentativo(i valori misurati che io do in input)
% tol: Tolleranza (precisione con cui voglio la soluzione)
% itmax: Numero massimo di iterazioni
%
% OUTPUT
% x: Soluzione
% iterations: Numero di iterazioni eseguite
% f: Vettore dei reidui (equazioni di chiusura)
% resid: Norma del vettore dei residui
% exitflag:  1 L'iterazione converge sul valore x, 0
%            0 Si e' superato il massimo numero di iterazioni

% Inizializza le variabili
resid=1000;
iterations=0;
while(resid>tol)
    % Calcola matrice Jacobiana
    J=feval(jacobian,x0);
    % Calcola residui
    f=feval(fun,x0);
    % Calcola norma del vettore dei residui
    resid=norm(f);
    
    % Trova il nuovo valore di x con il metodo di Newton-Raphson
    x = x0-J\f; % inv(J)  J \ 
    
    
    % Controllo sulla tolleranza
    if(resid<tol)
        exitflag=1;
        return
    end
    % Controllo sul numero delle iterazioni
    if(iterations>itmax)
        exitflag=0;
        return
    end
    % Assegan il nuovo valore come condizioni iniziali per la successiva
    % iterazione
    x0=x;
    iterations=iterations+1;
end

modulo main

clc; clear; close all;
x0=[257.85 ; 348.84]; 


str = sprintf('Soluzione di primo tentativo: %.1f  %.1f',x0);
disp(str);

% Criterio convergenza
tol=1e-16;
str = sprintf('Tolleranza: %.2e',tol);
disp(str);

% Massimo numero di iterazioni
itmax=100000;
str = sprintf('Massimo numero di iterazioni: %i',itmax);
disp(str);


str= sprintf ( '\n\n----Metodo Newton-Raphson----\n\n');
disp(str);
% Metodo di Newton-Raphson per trovare la soluzione del sistema di
% equazioni non lineari
[x,iterazioni,funzioni,residui,exitflag]=newtonraphson(@fun,@jacobian,x0,tol,itmax);


str = sprintf('Soluzione = %.2e %.2e',x);
disp(str);
str = sprintf('Iterazioni = %.i',iterazioni);
disp(str);
str = sprintf('Residui della funzione = %.3e %.3e',funzioni);
disp(str);
str = sprintf('Massima norma dei residui = %.3e',residui);
disp(str);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%                    PLOT                              %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r(2)=1;
r(3)=3;
r(4)=4;
r(6)=4;
r(7)=3;
%theta(2)=pi/4;
%theta(4)= pi;
for i=1:1:360   %qui come prima non sono sicuro sia giusto
    theta(1)=i
end
f(1)= r(2)*cosd( theta(1)) + r(3)*cosd( x(1)) + r(4)*cosd( x(2)) - r(6);
f(2)= r(2)*sind( theta(1)) + r(3)*sind( x(1)) + r(4)*sind( x(2)) + r(7);
z_point= r(2)*cosd( theta(1)) + r(3)*cosd( x(1)) + r(4)*cosd(x(2)) + r(6);




figure('Position',[100 100 900 500])
colormap jet % jet copper cool
title('\fontsize{15}Grafico delle posizioni')
grid on
xlabel('Theta1','fontweight','b');
ylabel('theta3 e theta2','fontweight','b');
hold on
plot( theta(1), x(1),'ro','Linewidth',4)
plot(theta(1), x(2), 'ro','Linewidth',4)
Top

[apatriarca]

Ad una rapida occhiata osservo che il codice

for i=1:1:360   %qui ho sicuramente sbagliato qualcosa perchè in output stampa i 360 valori ti theta
    theta(1)= i;
end

è equivalente a

theta(1) = 360;

se escludiamo la stampa dei vari valori di theta. Dubito sia quello che volevi ottenere.

Non mi è poi chiaro perché dici che ci sono 4 equazioni in 4 incognite.. Io vedo solo due equazioni ( \(f(1)\) e \(f(2)\) ) e 3 incognite (se si considera anche l'angolo \(\theta\)).

[andrea_d11]

Grazie per la risposta e per la pazienza che hai impiegato a guardare il codice.
Dunque il programma è parziale infatti qui trovo le prime due incognite noto che sia l'angolo theta (cioè theta varia tra 0 e 360, ed è questo che non riesco a far capire al programma, cioè di prendermi uno per volta i valori di \theta e per ognuno calcolarmi le due incognite), le altre due equazioni usano i risultati ottenuti da queste per calcolare le altre due incognite (che qui ho omesso per alleggerire la cosa).
Per quanto riguarda la tua osservazione sul ciclo for per definire theta, beh si era un tentativo disperato perchè in realtà non so come fare, avresti qualche suggerimento? Ci sono degli errori che non mi fanno arrivare al plot (scopo di tutto) finale con una serie di valori ma bensì con dei valori puntuali.