% Scripts for MAT 331 October 4, 2017 % Corrected after class % function to find zero by bisection function [x,k]=bisect(f,a,b) if sign(a)==sign(b) if sign(a)==0 x=a; k=0; return; else return; end k=0; while abs(b-a) > eps*abs(b) x = (a + b)/2; if sign(f(x)) == sign(f(b)) b = x; else a = x; end [a,b]; k = k + 1; end return; %--------------------------------------------------------- % function to apply Newton's method to a function % derivative must be provided function [y,k] = newton(f,fprime,x) k = 0; xprev=x+1; while abs(x - xprev) > eps*abs(x) xprev = x; x = x - f(x)/fprime(x); k = k + 1; end y=x; return; %--------------------------------------------- % function to apply secant method. Function and % two starting points must be given function [out,k] = secant(f,a,b) k=0; while abs(b-a) > eps*abs(b) c = a; a = b; b = b + (b - c)/(f(c)/f(b)-1); k = k + 1; end out=a; return; %--------------------------------------------------------- % function to apply Newton's method to a polynomial % derivative is computed by this function function [y,k] = poly_newton(p,x) p_prime=polyder(p); k = 0; xprev=x+1; while abs(x - xprev) > eps*abs(x) xprev = x; x = x - polyval(p,x)/polyval(p_prime,x); k = k + 1; if k>50 % stop if doesn't converge y=NaN; return; end % if end % while y=x; return; % try p=[1,0,0,1] cube roots of 1 % try p=[2,8,2,0,3,1] %-------------------------------------------------------- % function to plot result of Newton's method. Colors plane % according to which root that initial value finds % inputs % p = polynomial in MATLAB format: coefficeints in decreasing order % a = size box, picture will be of [-a,a]^2 % n = number of pixels along each edge of box function plot_newton(p,a,n) all_roots=roots(p); delta=abs(2*a)/n; % save_root=zeros(n); for k=1:n for j=1:n z=(-abs(a)+2*abs(a)*j/n)+1i*(-abs(a)+2*abs(a)*k/n); [y,iter]=poly_newton(p,z); if isnan(y) save_root(k,j)=0; else [m,root_num]=min(abs(all_roots-y)); save_root(k,j)=root_num; end end k end all_roots; save_root; figure; title('Domains of attraction for Newtons method', p) image(save_root,'CDataMapping','scaled') %axis image %colormap(flipud(jet(length(all_roots)))) toc; return; %-------------------------------------------- % function to count and find zeros function [counter,savezero]=count_zeros(f,a,b,n) delta=(b-a)/n; t=[a:delta:b]; counter=0; clear savezero; for k=1:n-1 if f(t(k))*f(t(k+1)) <= 0 counter=counter+1; savezero=fzero(f,[t(k),t(k+1)]); end end return;