Càlcul dels zeros d'una funció:

En aquesta pràctica veurem com es poden calcular numèricament els zeros d'una funció.

Teorema de Bolzano: Canvis de signe

Donada una funció f(x), per calcular els seus zeros haurem de fer:
% funció f(x) definida com a inline
f=@(x) sin(x);
%rang de valors que estudiem
xmin=-5;
xmax=5;
pas=0.15;
x=xmin:pas:xmax;
%valors de la funció
y=f(x);
 
%-------------------------------------------------------------
% pintem la funció per tenir una idea dels zeros
%-------------------------------------------------------------
semieix_y=1.1*max(abs(y)); %ajustem els eixos als valors de la funció
plot(x,y,'.-');
hold on;
axis([xmin xmax -semieix_y semieix_y]);
grid
line([xmin xmax],[0 0],'Color','k','LineWidth',2)
line([0 0],[-semieix_y semieix_y],'Color','k','LineWidth',2)

Determinació de l'interval de canvi de signe

%-------------------------------------------------------------
% Determinem els possibles zeros de la funció.
%-------------------------------------------------------------
%
% Fem un estudi del numero de zeros a partir dels canvis de signe de la
% funció
 
nzeros=0;
intervals=[];
for j=1:length(y)-1
if ( y(j)*y(j+1) <= 0 )
nzeros=nzeros+1;
plot(x(j),y(j),'ro')
plot(x(j+1),y(j+1),'yd')
disp(['zero ', num2str(nzeros)]);
disp(['interval [' num2str(x(j)) ',' num2str(x(j+1)) ']']);
intervals=[intervals; x(j), x(j+1)]; %guardem els intervals en una matriu
end
end
zero 1
interval [-3.2,-3.05]
zero 2
interval [-0.05,0.1]
zero 3
interval [3.1,3.25]

Mètode de la Bisecció

Es basa en dividir succesivament la meitat l'interval inicial [a,b], en el que hi ha un canvi de signe i anar descartant la meitat que ja no conté el canvi de signe (trobareu el codi al final del guió).
Iteracions de la bisecció

Càlcul del zero a cada interval

Ara per cada interval cridarem una funció que calculi el zero fins a una determinada tolerància d'error . Aquest criteri és suficient si la funció és prou vertical però en general és més precís utilitzar que que sempre donarà un valor més proper al zero real.
tol = 1.e-8; % demanem 10^(-8) de precisió en el càlcul
numZeros = size(intervals,1); %cada fila és un interval [xa, xb] de canvi de signe
for i = 1:numZeros
xa = intervals(i,1);
xb = intervals(i,2);
z(i) = biseccio(f, xa, xb, tol); %s'ha de triar el mètode a utilitzar
%z(i) = secant(f, xa, xb, tol);
%z(i) = newton(f,df,x0,tol);
end
iter = 1 x= -3.1250000000000000e+00 f(x)= -1.659189e-02 iter = 2 x= -3.1625000000000001e+00 f(x)= 2.090582e-02 iter = 3 x= -3.1437499999999998e+00 f(x)= 2.157345e-03 iter = 4 x= -3.1343749999999999e+00 f(x)= -7.217591e-03 iter = 5 x= -3.1390624999999996e+00 f(x)= -2.530151e-03 iter = 6 x= -3.1414062499999997e+00 f(x)= -1.864036e-04 iter = 7 x= -3.1425781250000000e+00 f(x)= 9.854713e-04 iter = 8 x= -3.1419921874999996e+00 f(x)= 3.995339e-04 iter = 9 x= -3.1416992187499995e+00 f(x)= 1.065652e-04 iter = 10 x= -3.1415527343749998e+00 f(x)= -3.991921e-05 iter = 11 x= -3.1416259765624996e+00 f(x)= 3.332297e-05 iter = 12 x= -3.1415893554687497e+00 f(x)= -3.298121e-06 iter = 13 x= -3.1416076660156245e+00 f(x)= 1.501243e-05 iter = 14 x= -3.1415985107421873e+00 f(x)= 5.857152e-06 iter = 15 x= -3.1415939331054688e+00 f(x)= 1.279516e-06 iter = 16 x= -3.1415916442871090e+00 f(x)= -1.009303e-06 iter = 17 x= -3.1415927886962889e+00 f(x)= 1.351065e-07 iter = 18 x= -3.1415922164916990e+00 f(x)= -4.370981e-07 iter = 19 x= -3.1415925025939941e+00 f(x)= -1.509958e-07 iter = 20 x= -3.1415926456451415e+00 f(x)= -7.944652e-09 -------------------------------- iter = 1 x= 2.5000000000000355e-02 f(x)= 2.499740e-02 iter = 2 x= -1.2499999999999734e-02 f(x)= -1.249967e-02 iter = 3 x= 6.2500000000003109e-03 f(x)= 6.249959e-03 iter = 4 x= -3.1249999999997113e-03 f(x)= -3.124995e-03 iter = 5 x= 1.5625000000002998e-03 f(x)= 1.562499e-03 iter = 6 x= -7.8124999999970579e-04 f(x)= -7.812499e-04 iter = 7 x= 3.9062500000029698e-04 f(x)= 3.906250e-04 iter = 8 x= -1.9531249999970440e-04 f(x)= -1.953125e-04 iter = 9 x= 9.7656250000296291e-05 f(x)= 9.765625e-05 iter = 10 x= -4.8828124999704056e-05 f(x)= -4.882812e-05 iter = 11 x= 2.4414062500296117e-05 f(x)= 2.441406e-05 iter = 12 x= -1.2207031249703969e-05 f(x)= -1.220703e-05 iter = 13 x= 6.1035156252960739e-06 f(x)= 6.103516e-06 iter = 14 x= -3.0517578122039478e-06 f(x)= -3.051758e-06 iter = 15 x= 1.5258789065460631e-06 f(x)= 1.525879e-06 iter = 16 x= -7.6293945282894233e-07 f(x)= -7.629395e-07 iter = 17 x= 3.8146972685856038e-07 f(x)= 3.814697e-07 iter = 18 x= -1.9073486298519098e-07 f(x)= -1.907349e-07 iter = 19 x= 9.5367431936684699e-08 f(x)= 9.536743e-08 iter = 20 x= -4.7683715524253140e-08 f(x)= -4.768372e-08 iter = 21 x= 2.3841858206215780e-08 f(x)= 2.384186e-08 iter = 22 x= -1.1920928659018680e-08 f(x)= -1.192093e-08 iter = 23 x= 5.9604647735985499e-09 f(x)= 5.960465e-09 -------------------------------- iter = 1 x= 3.1750000000000007e+00 f(x)= -3.340113e-02 iter = 2 x= 3.1375000000000006e+00 f(x)= 4.092642e-03 iter = 3 x= 3.1562500000000009e+00 f(x)= -1.465682e-02 iter = 4 x= 3.1468750000000005e+00 f(x)= -5.282322e-03 iter = 5 x= 3.1421875000000004e+00 f(x)= -5.948464e-04 iter = 6 x= 3.1398437500000007e+00 f(x)= 1.748903e-03 iter = 7 x= 3.1410156250000005e+00 f(x)= 5.770286e-04 iter = 8 x= 3.1416015625000004e+00 f(x)= -8.908910e-06 iter = 9 x= 3.1413085937500007e+00 f(x)= 2.840598e-04 iter = 10 x= 3.1414550781250004e+00 f(x)= 1.375755e-04 iter = 11 x= 3.1415283203125002e+00 f(x)= 6.433328e-05 iter = 12 x= 3.1415649414062505e+00 f(x)= 2.771218e-05 iter = 13 x= 3.1415832519531257e+00 f(x)= 9.401637e-06 iter = 14 x= 3.1415924072265629e+00 f(x)= 2.463632e-07 iter = 15 x= 3.1415969848632814e+00 f(x)= -4.331273e-06 iter = 16 x= 3.1415946960449221e+00 f(x)= -2.042455e-06 iter = 17 x= 3.1415935516357427e+00 f(x)= -8.980459e-07 iter = 18 x= 3.1415929794311528e+00 f(x)= -3.258414e-07 iter = 19 x= 3.1415926933288576e+00 f(x)= -3.973906e-08 iter = 20 x= 3.1415925502777102e+00 f(x)= 1.033121e-07 iter = 21 x= 3.1415926218032837e+00 f(x)= 3.178651e-08 iter = 22 x= 3.1415926575660706e+00 f(x)= -3.976277e-09 --------------------------------
disp('zeros :');
zeros :
z
z = 1×3
-3.141592645645142e+00 5.960464773598550e-09 3.141592657566071e+00

Exercici 1:

Canvieu el criteri de detecció del zero de bisecció per la condició .
Hint: En aquest cas només cal fer i només calen un parell més d'iteracions.
Agafeu la funció en que és una funció molt plana (tota ella està per sota de ) i mireu la diferència entre els dos criteris (el segon sempre és més precís!!)

Exercici 2: Mètode de la secant (regula falsi)

Implementeu el mètode de la secant creant una funció de Matlab
function zero = secant(f, xa, xb, tol)
i substituiu en el programa anterior la crida a biseccio per la crida a secant.
Mètode de la secant:
Anàlogament al que es fa en el mètode de bisecció, es pot fer el mateix però triant com a punt intermedi c, el punt de tall de l'eix x (y=0) amb la recta que uneix els dos punts (a,f(a)) i (b,fb)).
Podeu comprovar que la fórmula és:
A partir d'aquí es procedeix exactament igual que a bisecció, es descarta el subinterval que no conté el zero.
(source https://planetcalc.com/3712/)
Solució: Pel primer dels zeros només calen dues iteracions i donen:
iter = 1 x= -3.1415994094031583e+00 f(x)= 6.755813e-06
iter = 2 x= -3.1415926532799121e+00 f(x)= -3.098811e-10

Exercici 3: Mètode de Newton

Implementeu el mètode de Newton creant una funció de Matlab
function zero = newton(f, df, x0, tol)
i substituiu en el programa anterior la crida a biseccio per la crida a newton.
Mètode de Newton:
Anàlogament al que es fa en el mètode de secant, es pot fer el mateix però triant com a punt intermedi c, el punt de tall de l'eix x (y=0) amb la recta tangent en un punt inicial de l'interval [a,b]. Per exemple el punt mig de [a,b], .
Això requereix, entre d'altres coses, el càlcul de la derivada de la funció f(x) i definir-la també com una funció inline.
La fórmula matemàtica del mètode en aquest cas és:
A partir d'aquí es procedeix exactament igual que a bisecció, es descarta el subinterval que no conté el zero.