Raíces de una ecuación

La ecuación de segundo grado

a x 2 +bx+c=0⇒{ x 1 = -b+ b 2 -4ac 2a x 2 = -b- b 2 -4ac 2a

Creamos el script

a=input('primer coeficiente');
b=input('segundo coeficiente');
c=input('tercer coeficiente');
x1=(-b+sqrt(b^2-4*a*c))/(2*a);
x2=(-b-sqrt(b^2-4*a*c))/(2*a);
disp([x1 x2])
disp(x1)
disp(x2)

Lo guardamos con el nombre ecuacion

También podemos crear una función:

function [x1,x2] = ecuacion_segundo(a,b,c)
    x1=(-b+sqrt(b^2-4*a*c))/(2*a);
    x2=(-b-sqrt(b^2-4*a*c))/(2*a);
 
end

Lo guardamos con el nombre ecuacion_segundo. Repitiendo los ejercicios anteriores:

[z1 z2]=ecuacion_segundo(1,4,6)
    z1 = -2.0000 + 1.4142i
    z2 = -2.0000 - 1.4142i

>> [v1 v2]=ecuacion_segundo(1,-2,-3)
    v1 = 3
    v2 = -1

Ecuaciones de grado 'n'

Se emplea la función roots, se pone un vector con los coeficientes de la ecuación y nos da las soluciones (en columna).

>> p=[an, an-1... a1, a0]
>> x=roots(p)

NOTA: con i y -i se aconseja poner 1.

>> a=[1i -1i];
>> poly(a)
    ans =     1     0     1

La función polyval permite calcular el valor numérico de un polinomio para un determinado valor de la variable x.

Ejemplo: f(x)=x5+5x4-3x3-15x2-4x-20 , si queremos hallar f(11)

>> p=[1 5 -3 -15 -4 -20];
>> polyval(p,11)
    ans =  228384
>> polyval(p,0)
    ans =    -20

Esto nos permitirá conocer si ciertos valores son raíces o no de la ecuación. Comprobando con el ejercicio anterior:

>> polyval(p,-5)
    ans =   0
>> polyval(p,i)
    ans =   0
>> polyval(p,-i)
    ans =   0
>> polyval(p,-2)
    ans =  0
>> polyval(p,2)
    ans =   0

Raíces, utilizando solve de Math Symbolic

>> syms x;
solve(x^3+6*x^2+11*x+6,x) 
    ans = -3,-2,-1

Comprobamos con roots

>> p=[1 6 11 6];
>> roots(p)
    ans =   -3.000
            -2.0000
            -1.0000

Ejemplo:

>> a=[1 1 1];
>> roots(a)
    ans =  -0.5000 + 0.8660i
            -0.5000 - 0.8660i
>>  solve(x^2+x+1,x)
    ans =  - 1/2 + (3^(1/2)*i)/2
        - 1/2 - (3^(1/2)*i)/2
>> double(ans)
    ans =  -0.5000 - 0.8660i
            -0.5000 + 0.8660i

En vez de double(ans) se puede poner vpa(ans) o vpa(ans,k) donde k es el número de cifras

Ejercicios

  1. Hallar las raíces de las siguientes ecuaciones:

    2. { 3 x 2 -2x-8=0 x 4 -8 x 3 +26 x 2 -40x+25=0 x 4 +4 x 3 +10 x 2 +12x+5=0  { x 5 +4 x 4 +5 x 3 -6x-4=0 x 4 +13 x 2 +36=0 x 3 +6 x 2 +11x+6=0

  2. En las dos primeras, comprobar que los valores obtenidos son las soluciones de las ecuaciones.

  3. Si la gráfica de función y=f(x) corta al eje OX en los puntos x=-4,1,3 y 7. Calcular f(21).

Ecuaciones transcendentes

En la ingeniería a menudo nos encontramos con ecuaciones que no se pueden resolver por los métodos analíticos, para estos casos intentaremos obtener una solución aproximada.

El teorema de Bolzano, que nos permitirá localizar las soluciones afirma que si una función y=f(x) es continua en un intervalo cerrado [a, b] y toma valores de signo opuesto en los extremos del mismo, dicha función se anula, por lo menos en un punto interior.

En este apartado trabajaremos con funciones anónimas que nos permiten definir una función sin necesidad de guardarla en un fichero .m. Ejemplos

>> f=@(x) 3*x-exp(x);
>> g=@(x) 5*x^2-exp(x);
>> h=@(x) x-cos(x);
>> j=@(x) x/3-sin(x);
>> k=@(x) x-exp(sin(x));

En el Workspace aparecen con el siguiente símbolo

Método gráfico

Tenemos la ecuación : 3x-ex=0.

y=f(x)=3x-ex es continua en todo R

>> f(0)
    ans =    -1
>> f(1)
    ans =    0.2817
>> f(2)
    ans =   -1.3891

Con fplot

>> f=@(x)3*x-exp(x);
>> fplot(f,[0 2])
>> grid on

Nos acercamos a la primera raíz

>> [x y]=ginput
    x =    0.6191
    y =  -1.2062e-11

Con la otra raíz:

>> [x y]=ginput
    x =    1.5121
    y =  -5.3180e-13

De otra forma: Intersección

>> n=@(x)exp(x);
>> fplot(n,[0 2])
>> grid on
>> hold on
>> m=@(x) 3*x;
>> fplot(m,[0 2])
>> hold off

Nos acercamos a la primera intersección

>> [x,y]=ginput
    x =    0.6191
    y =    1.8572

>> [x,y]=ginput
    x =    1.5121
    y =    4.5364

Ahora utilizando la función fzero de Matlab:

>> fzero(f,[0 1])
    ans =    0.6191
>> fzero(f,[1 2])
    ans =    1.5121

Nota.- La función fzero tiene que cambiar de signo en el intervalo seleccionado, caso contrario se produce un error.

Ejemplos

Método de las aproximaciones sucesivas

Supongamos la ecuación f(x)=0, siendo f(x) una función continua en el dominio donde estamos trabajando.

Sustituimos f(x)=0 por una ecuación equivalente x=g(x) (no siempre se puede hacer). En el ejemplo anterior x= e x 3 ⇒g(x)= e x 3

Suponiendo un valor aproximado de la raíz x0, lo sustituimos en g(x) y obtenemos x1=g(x0), lo volvemos a sustituir y tenemos x2=g(x1) y así, sucesivamente, xn=g(xn-1)

Si existe lim   x n n→∞ =p , p es la raíz buscada

Elaboramos el siguiente script

x0=input('valor inicial: ');
iter=input('número de iteraciones: ');
for i=1:iter
    x=exp(x0)/3;
    x0=x;
end
disp(x)

y lo guardamanos en el fichero trascendente, tenemos:

>> trascendente
valor inicial: 0
número de iteraciones: 12
    0.6180

>> trascendente
valor inicial: 0
número de iteraciones: 20
    0.6190

Vemos que depende del número de iteraciones que le ponemos.

Nota.-También podemos poner x=f(x0)) (genérico) si previamente hemos puesto

f=@(x) exp(x)/3;
x0=input('valor inicial: ');
iter=input('número de iteraciones: ');
for i=1:iter
    x=f(x0);
    x0=x;
end
disp(x)

Sin embargo, si en vez de utilizar un script genérico lo hacemos como una función genérica, hay que hacer llamada a la función ejemplo: x=transc(f,x0,a).

function  x  = transc(f,x0,a)
    for i=1:a
        x=f(x0);
        x0=x;
        end
end

Sabemos que tiene otra solución x= 1.5121. Por mucho que pongamos el valor inicial cercano

>> trascendente
valor inicial: 1.5
número de iteraciones: 50
    0.6191

Con este método, si hemos encontrado g(x) a partir de la ecuación que nos daban, pero cuando la derivada de g(x) en valor absoluto en un entorno de la solución es mayor que la unidad no nos sirve, pues la sucesión {xn} diverge (no sirve si es >1).

g(x)= e x /3⇒g'(x)= e x /3

Poniendo h(x)=g'(x):

>> h=@(x) exp(x)/3;
>> h(1.4)
    ans =    1.3517

No lo cumple, sin embargo

>> h(0.5)
    ans =    0.5496

Si tomamos otro ejemplo: x-cosx=0.

g(x)=cosx, g'(x)=-senx, cumple la condición. Tenía una solución 0.7391

>> g=@(x) cos(x);
>> x=transc(g,0,10)
    x =    0.7314
>> x=transc(g,0,20)
    x =    0.7389
>> x=transc(g,0,40)
    x =    0.7391

Vemos las diferencias según el número de iteraciones que ponemos. Debido a esto, tendremos que poner una condición de FINALIZACION para obtener la solución con una precisión determinada. Mediante el error relativo:

| x n - x n-1 x n |<ε  

Para ello sustituimos el bucle for por un bucle while indefinido, que terminará cuando se cumpla la condición.

x0=input('valor inicial: ');
Error=0.00001;
while(1)
    x=exp(x0)/3;
    if abs((x-x0)/x)<Error
        break
    end
        x0=x;
end
disp(x)

Nota: while(1) siempre se ejecuta

Lo guardamos en el fichero trascendente1. Tenemos

>> trascendente1
valor inicial: 0
    0.6191

También lo podemos hacer con una función (poniendo diferentes valores iniciales y errores).

function x = raiz(x0, a)
    while(1)
        x=exp(x0)/3;
        if abs((x-x0)/x)<
            break
        end
        x0=x;
    end
end

La llamada

>> raiz(0,0.00001)
    ans =    0.6191

Alternativamente, con un script

x0=input('valor inicial: ');
Error=0.00001;
while(1)
    x=exp(x0)/3;
    if abs((x-x0)/x)<Error
        break
    end
    x0=x;
end
disp(x)

Lo guardamos en el fichero raiz2.

>> raiz2
valor inicial: 0
    0.6191

Podemos generalizar el método (no dependa de g(x)) de la siguiente manera:

function x=generalizacion(f,x0,error)
    while (1)
        x=f(x0);
        if abs((x-x0)/x)<error
            break
        end
        x0=x;
    end
end

En la ventana de Command Window definimos f y lo aplicamos

>> f=@(x) exp(x)/3;
>> generalizacion(f,0,0.000001)
    ans =    0.6191
>> g=@(x) cos(x);
>> generalizacion(g,0,0.000001)
    ans =    0.7391

Ejercicios

  1. x/4-lnx=0

  2. 2x-e-senx=0

  3. x-5ln(3x)=0

Soluciones

Método de Newton-Raphson

Si desarrollamos la función f(x+h) en serie de Taylor en el punto x:

f(x+h) =f(x)+f'(x)·h+...

Si h es pequeño podemos despreciar los términos en, h2, h3... Si x+h es raíz de la ecuación

f(x+h)=0⇒h≈- f(x) f'(x)

Creamos un desarrollo iterativo que nos permita calcular la raíz de la siguiente forma:

x n+1 = x n +h= x n - f( x n ) f'( x n )

Definimos el procedimiento de Newton-Raphson de forma análoga al anterior de las aproximaciones sucesivas.

function x=newton_raphson(f,f_prima,x0,error)
    while(1)
        x=x0-f(x0)/f_prima(x0);
        if abs((x-x0)/x)<error
            break
        end
        x0=x;
    end 
end

Repitiendo el mismo ejemplo:

>> f=@(x) x-exp(x)/3;
>> f_prima=@(x) 1-exp(x)/3;
>> format long
>> newton_raphson(f,f_prima,0.5,0.00001)
    ans =   0.619061286735945

Para la otra raíz:

>> newton_raphson(f,f_prima,1.4,0.000001)
    ans =   1.512134551659757

Ejemplos

Haciendo con un bucle for

function x = Rahfor(f,f_prima,x0,iter)
    for i=1:iter
        x=x0-f(x0)/f_prima(x0);
        x0=x; 
    end
end
>> x=Rahfor(f,f_prima,1,100)
    x =   1.429611824725556
>> x=Rahfor(f,f_prima,8,100)
    x =   8.613169456441398

Método del punto medio o bisección

Este método se basa en el teorema de Bolzano que nos dice:

Si una función f(x) es continua en un intervalo [a, b] y en los extremos del intervalo toma valores de signo opuesto, entonces existe, al menos, un valor c interior al intervalo donde f(c)=0.

Tenemos la ecuación f(x)=0, donde la función f(x) cumple las condiciones de Bolzano.

m= a+b 2 ⇒{ f(m)=0⇒m='c'  raiz buscada. f(a)·f(m)<0  'c',estará en el intervalo (a, m). f(a)·f(m)>0  'c', estará en el intervalo (m, b).

Repitiendo el proceso obtenemos una sucesión de intervalos:

[ a 1 , b 1 ]⊃[ a 2 , b 2 ]⊃....[ a n , b n ]⊃... tal que: f( a n )·f( b n )<0 y la amplitud  b n − a n = 1 2 n (b-a) lim n→∞    a n = lim n→∞    b n =c

Nos podemos plantear el error de cálculo de dicha raíz. Si queremos que sea menor que ε

2 n > b-a ε ⇒n> 1 ln2 ·ln( b-a ε ).

Condiciones de terminación del proceso

function m = mediofor(f,a,b,error,iter)
    for i=1:iter
        m=(a+b)/2;
        if abs((b-a)/m)<error
            break
        elseif f(a)*f(m)<0
            b=m;
        else
            a=m;
        end
    end
end

Aplicado a nuestro ejemplo:

>> f=@(x) 3*x-exp(x);
>> m= mediofor(f,-1,1,1e-10,20)
    m =   0.619062423706055
>> m= mediofor(f,-1,1,1e-10,200)
    m =   0.619061286764918
>> fzero(f,[-1,1])
    ans =   0.619061286735945

Para el otro punto

>> m=mediofor(f,1,2,1e-10,200)
    m =   1.512134551710915
>> m=mediofor(f,1,2,1e-10,1000)
    m =   1.512134551710915
>> fzero(f,[1,2])
    ans =   1.512134551657843

Con un bucle while

function m = mediowh( f,a,b,error )
    while(1)
        m=(a+b)/2;
        if abs((b-a)/m)<error
             break
        elseif f(a)*f(m)<0
          b=m;
        else
          a=m;
        end
    end  
end
>> m=mediowh(f,-1,1,1e-20)
    m =   0.619061286735946		0.619061286735945(fzero)
>> m=mediowh(f,1,2,1e-20)
    m =   1.512134551657844		1.512134551657843(fzero)

La diferencia es que con el while si la función en a y b no tiene signos contrarios, Matlab se cuelga. Con el for avisa que no ha encontrado... Puede ser porque el número de iteraciones no son suficientes o porque se ha afinado mucho con el error.

Para el otro ejercicio:

>> f=@(x) x/2-3*log(x);
>> m=punto_medio(f,1,2,40)
    m =   1.226867675781250
>> m=punto_medio(f,10,20,40)
    m =  16.998901367187500

Ejemplos

Newton-Raphson

Punto medio

fzero

>> fzero(u,[1 2])
    ans =   1.784791163424146
>> fzero(u,[3 4])
    ans =   3.634835618651561