%bsp_46
function x = trisect_newton(func,a,b,tol,maxiter);
% Doku
x = trisect(func,a,b,tol);
[x,X] = newton(func,x,maxiter);
if length(X) >=maxiter
    tol = input('decrease tol for trisection, tol = ' )
    x = trisect_newton(func,a,b,tol,maxiter);
end

% bsp_47
function [f,df] = func(x);
%
f = sign(x-1).*abs((x.^11-1)).^(1/11) + .5 + .05*sin(x/100);
df = (abs(x.^11-1)).^(-10/11) .*x.^(10) + 5e-4*cos(x/100);


%bsp_48
function [r,c] = rate(x);
%
% Konvergenzrate von x
%
diff = log(abs(x(2:end)-x(1:end-1)));

p=polyfit(diff(1:end-1),diff(2:end),1);
%plot(diff(1:end-1),diff(2:end),'*')
r = p(1);
c = exp(p(2));
return


%bsp_50
n=1000;c = 1.13198824;
x = zeros(n,1);
x(1) = 1; x(2) = 2;
for i = 3:n
    x(i) = x(i-1) + 2*(round(rand)- .5) * x(i-2);
end
plot(log(abs(x)));
hold on
xn = 1:n;
plot( xn, xn*log(c),':')

%bsp_51
function x1 = runde(z,n);
% rounds elements of inputarray to n decimal digits 
% use: x= runde(z,n);
%

[row,col] = size(z);

ind = find(z);
y = z(ind);

signy= sign(y);
absy= abs(y);
logy= log10(absy); 

% normalized reprsentation: 0,abc... 10^k, a>0
logy = ceil(logy);

absy = absy ./ 10.^(logy);

% round
absy= absy .* 10^n;
x = round(absy)/ 10^n;
x= signy .* x .* 10.^logy;
x1 = zeros(size(z));
x1(ind) = x;

return

%bsp_52
function sum =add(x,y,n)
% simulation of addition with n decimal digit accuracy
% needs function runde
% length of mantissa of input n!!

sum = runde(x+y,n) ;


% bsp 53
function x=backsubstitute(a,b);
% backsubstitution for upper triangular system
%
if any(abs(diag(a)) <= 1e-13)
    error('matrix close to singular')
end 
n=length(b);

x=zeros(n,1);
x(n)=b(n)/a(n,n);
for i=n-1:-1:1
   x(i) = (b(i) - a(i,i+1:n)*x(i+1:n))/a(i,i);
end

% Bsp_54
function [a,l]=gauss(a);
% Gauss elimination without interchanges

tol = 1e-14;
[row,col] = size(a);
l=eye(row);

for k=1:min(row,col)-1
   if(abs(a(k,k))<tol)
      disp('elimination not possible')
      return
   end
   
   z=a(k,k+1:col);
   
 %  for i=k+1:row
 %     m=a(i,k)/a(k,k);
 %     a(i,k+1:col) = a(i,k+1:col)-z*m;
 %  end
 m=a(k+1:row,k)/a(k,k);
 a(k+1:row,k+1:col) = a(k+1:row,k+1:col) -m*z;
 a(k+1:row,k)= zeros(row-k,1);%
% l(k+1:row,k)=m;
end


% bsp_55
 A=[0,0,0]; B=[1,-1,3]; C=[-1,1,1];
a = sqrt((B-C)*(B-C)');
b = sqrt((C-A)*(C-A)');
c = sqrt((B-A)*(B-A)');

disp(['Die Länge der Seite a beträgt ',num2str(a),' cm'])
fprintf('Die Länge der Seite b beträgt %5.2f cm\n', b)
fprintf('Die Länge der Seite c beträgt %-5.2f cm\n', c)



clear

% bsp_56

x= .5:.1:1; y= log10(x);
fid = fopen('Logarithmen','w');
fprintf(fid,' x   log(x)\n\n');
fprintf(fid,'%3.1f  %7.4f\n',[x;y]);
fclose(fid);

x=1.1:.1:1.5;y=log10(x);
fid = fopen('Logarithmen','a');
fprintf(fid,'%3.1f  %7.4f\n',[x;y]);
fclose(fid);

% Einlesen der Daten
fid = fopen('Logarithmen','r');
s=fscanf(fid,'%s',2);    %überlese 1.Zeile, Leerzeile wird ignoriert
%fscanf(fid,'%s');
X = fscanf(fid,'%g',[2,inf]);
fclose(fid);
z= X(1,:)',y=X(2,:)'

% Bsp_59
clear
% fid =fopen('polygon_dat','r');
% X = fscanf(fid, '%g',[2,inf]);
% fclose(fid);
%
% alternativ
%
X = load('polygon_dat');
x= X(1,:);y=X(2,:);
fill(x,y,'y');
axis equal
hold on
plot([x,x(1)],[y,y(1)],'b')
hold off



clear 

% Bsp_58

x= newton2(@func,pi/2,50);
X = load('newton_dat');
[r,c] = rate(X(1:end-1,2));     % letzte Differenz = 0!
fprintf('Das Newtonverfahren konvergiert mit Rate r = %5.3f\n',r);


