%FY1001/TFY4145 Fysikk høsten 2014.
%Øving 3, oppgave 5: Numerisk bestemmelse av
%vinkler for fire like masser som henger med jevne
%mellomrom på ei snor.
%
clear all;
format long;
%N = max antall iterasjoner
N = 100;
%L = snoras lengde
L = 1.0;
%D = lengden mellom festepunktene
D = 0.9;
%Vi har utledet en ligning for x = cos(alpha):
%(L/5)*(1+2x+4x/sqrt(1+3x^2))=D.
%Denne ligningen kan omformes til en 4.gradsligning i x,
%og på sett og vis løses analytisk (se Rottmann). Men løsningen
%blir stygg og gir lite innsikt. Derfor vil vi løse problemet numerisk.
%En mulighet (det finnes flere!) er da å omforme ligningen til
%formen x = f(x), sette inn en startverdi
%x0 på høyre side, regne ut f(x0) og sette x1 = f(x0), deretter
%regne ut f(x1) og sette x2 = f(x1) osv. Hvis dette "iterative
%skjemaet" konvergerer, vil vi etter N iterasjoner med god tilnærmelse
%ha xN = x(N-1), som er løsningen. (Se Fixed-point iteration, wikipedia.)
%Her er flere varianter aktuelle. Vi innfører gamma=(5D/L-1)/2 og har:
%1: x = gamma/(1+2/sqrt(1+3x^2))
%2: x = gamma - 2x/sqrt(1+3x^2)
%3: x = sqrt(1+3x^2)(gamma-x)/2
gamma = (5*D/L - 1)/2;
%Setter startverdi for x: x0 = 0.5;
x0 = 0.5;
%Setter startverdi for feilen lik 1.0 slik at while-løkken garantert starter
error = 1.0;
%Setter max tillatt feil i x lik 1e-6
maxerror = 1e-6;
j = 0;
while (error > maxerror && j < N)
    x1 = gamma/(1+2/sqrt(1+3*x0*x0));  %Metode 1
    %x1 = gamma - 2*x0/sqrt(1+3*x0*x0); %Metode 2
    %x1 = (gamma-x0)*sqrt(1+3*x0*x0)/2; %Metode 3
    error = abs(x1-x0);
    x0 = x1;
    j = j+1;
end;
%Gi  beskjed dersom vi har brukt opp max antall iterasjoner
%uten at feilen har blitt mindre enn maxerror
if j == N;
    fprintf('Max antall iterasjoner brukt opp.\n');
else
    fprintf('Vinkelen alpha er %3.12f',acos(x0)*180/pi);
    fprintf(' grader.\n');
end;
alpha = acos(x0);
%Vi kan nå bestemme vinkelen beta, samt de 3 snordragene i enheter av mg:
beta = atan(tan(alpha)/2);
S1 = 2/sqrt(1-x0*x0);
S2 = sqrt((1+3*x0*x0)/(1-x0*x0));
S3 = 2*x0/sqrt(1-x0*x0);
fprintf('Vinkelen beta er %3.12f',beta*180/pi);
fprintf(' grader.\n');
fprintf('S1 = %3.12f\n',S1);
fprintf('S2 = %3.12f\n',S2);
fprintf('S3 = %3.12f\n',S3);
%Plotter til slutt funksjonen x og de tre
%variantene av høyresider f(x) mellom 0 og 1
x = 0:0.001:1;
hold on;
plot(x,x,'k');
plot(x,(gamma-x).*sqrt(1+3.*x.*x)./2,'b');
plot(x,gamma - 2.*x./sqrt(1+3.*x.*x),'g');
plot(x,gamma./(1+2./sqrt(1+3.*x.*x)),'r');
