Code for using the newton iteration method

Sorry about the french in the code, my class is in french. Starting with 2 equations, I isolated phi in relation to s and plugged it into my other equation. I have an equation for s that is equal to zero and an equation for phi which works fine. My problem is that when I use the newton iteration method to find the value of s, I cant find non imaginary values. The output keeps giving me imaginary numbers for s, and a lot of them are in the 200+ range and they shouldnt get bigger than like 20.

The function itself works fine for normal values but I keep getting imaginary s numbers which mess up my equation. This makes me think theres some kind of error in the way i define deltas, which is the variation of s.

I know I need to use the absolute value of deltas in the while but when I put it in it breaks the code and runs forever giving the same value

I also know the fzero function exists but I'm not supposed to use it.

SN = 12;

%une gamme qui fait chaque 5 degrés
theta=linspace(0,2*pi, 72);
phi=zeros(1, 72);
s=zeros(1, 72);
deltas=0;


%fonction qu'on trouve en remplacant tout les phis par des s
syms svar thetavar
f = 10*sin(acos((3*cos(thetavar)-svar)/10)) + 1/3*cos(svar/SN)-2/3*cos(svar*sqrt(SN)/10)^2-3*sin(thetavar);

%ds=(f(s)/(df/ds))
ds=f/diff(f, svar);


lim=0.001;
j=1;

%puisque phi, s et theta sont des array il faut faire un loop pour les
%mesurer
while(j<=72)

    j
    firstloop=1;

    while((deltas) > lim | firstloop==1)
        %variation de s (deltas) en insérant les valeurs courantes de s et
        %theta et il faut aussi le convertir en double

        deltas=double(subs(ds, {svar, thetavar}, {s(j), theta(j)}))

        %correction de s
        s(j)=s(j)+deltas;
        
        firstloop=0;
    end

    phi(j)=f2(theta(j), s(j));

    %démontre les résultats du code
    distance=s(j)
    j=j+1;

end

figure(1);
plot(theta, s);

title('s in relation to Theta')
xlabel('Angle (RAD)');
ylabel('Distance (cm)')
axis([0 2*pi -20 20])

figure(2);
plot(theta, phi);

title('Phi in relation to Theta')
xlabel('Angle (RAD)');
ylabel('Angle 2 (RAD)')
axis([0 2*pi 0 2*pi])


%l'équation qu'on insère s pour trouver phi
function f = f2(theta, s)
    f = acos(acos((3*cos(theta)-s)/10));
end

[Torsten]

The "acos(x)" function only gives real values for arguments x that are less or equal 1 in absolute value. I guess your iteration gives bigger values.

Answer 1

[Aman]

Hi Jean,

I understand that you are getting imaginary values instead of real numbers while updating "s" with the help of the Newton-Raphson method.

The issue that you are experiencing is that after the first iteration, the value of "deltas" is quite high, due to which the below part of your function gives out an imaginary value as "acos" receives value beyond the [-1, 1] range.

10*sin(acos((3*cos(thetavar)-svar)/10))

In order to rectify the issue, either correct your function or put a check so that "acos" receives a value in the range [-1, 1].

I hope this helps!