Inicio
MATLAB
Raíces de ecuaciones Sistemas de ecuaciones Valores y vectores propios Integración numérica Ecuaciones diferenciales Interpolación, Interpolación, regresión
Numérico
Funciones MATLAB para calcular las raíces de una ecuación En esta página, vamos calcular la raíces de un polinomio mediante la función roots roots.. Un estudio más detallado de los polinomios se encuentra en la página titulada "Polinomios " Polinomios y fracciones polinómicas". polinómicas ". La raíz de una ecuación mediante procedimientos gráficos, que nos va a servir para estimar un valor próximo a la raíz buscada. Las raíces de una ecuación transcendente mediante fzero mediante fzero y y finalmente, las raíces de un sistema de ecuaciones transcedentes mediante fsolve mediante fsolve..
Raíces de un polinomio Para calcular las raíces de la ecuación
Ejercicios
a1 xn+a2 xn-1+...+an x+an+1=0 se emplea la función roots roots y y se le pasa el vector p vector p formado formado por los coeficientes del polinomio. La función roots devuelve un vector columna que contiene las raíces. >> p=[a1 a2 ... an >> x=roots(p)
an+1];
La función roots tiene una función inversa poly inversa poly que que se le pasa el vector x vector x que que contiene las raíces y devuelve los roots tiene coeficientes del polinomio p=poly(( x p=poly x)) Sea f Sea f ( x)= x)= x x5- 3. 3 .5 x4 + 2.75 x3 + 2. 2.125 x2 - 3. 3.875 x + 1.25 Guardamos los coeficientes del polinomio en el vector fila p fila p=[1 =[1 -3.5 2.75 2.125 -3.875 1.25]; Mediante la función plolyval función plolyval , podemos calcular el valor del polinomio cuando proporcionamos el valor de x >> p=[1 -3.5 2.75 2.125 -3.875 1.25]; >> polyval(p, 1.5) %valor del polinomio cuando se proporciona el valor de x. ans= -0.6250 >> r=roots(p) %raíces del polinomio r = 2.0000 -1.0000 1.0000 + 0.5000i 1.0000 - 0.5000i 0.5000 >> c=poly(r) %reconstruimos el polinomio a partir de las raíces
Comprobamos que las raíces calculadas son correctas utilizando la función polyval, polyval, pasándole pasándole los coeficientes del polinomio p y polinomio p y el valor x valor xi de cada una de las raíces >> polyval([1 3 3 2],
-2)
Como ejercicio se sugiere al lector, hallar las raíces de las ecuaciones siguientes con la función roots roots y y reconstruir el polinomio con la función poly función poly,, o calcular el valor del polinomio para cada una de las raíces con la función polyval x3 + + x x2 + + x x + 1 = 0 5 2 x +x -2 x x-4=0 -4=0 5 4 3 x - x x + + x x +2 x2-1
Procedimiento gráfico Vamos a estudiar varios procedimientos para calcular la raíz de una ecuación trascedente, por ejemplo, xcos( x)=0 En primer lugar, vamos a hacer una representación gráfica de esta función para conocer aproximadamente el punto donde la función corta al eje X. Escribimos el script mouse _ raiz y lo guardamos un fichero .M x=0:0.1:pi/2; y=x-cos(x); plot(x,y) grid on xlabel('x') ylabel('y') title('x-cos(x)')
Se selecciona en el menú Tools/Data Cursor, o el icono de la barra horizontal de herramientas . Situamos el cursor en forma de cruz próximo a la raíz de la función y podemos leer sus coordenadas. Cuanto mejor sea la resolución de la gráfica más cerca podremos estar de la raíz buscada. La resolución de la gráfica es el paso Δ x utilizado, para calcular la tabla de valores ( x, y) en el intervalo que va desde xi hasta x f . En este ejemplo Δ x=0.1
Podemos utilizar la herramienta Zoom In y Pan para acercarnos a este punto y el botón Zoom Out para alejarnos y volver a la situación inicial. En la imagen podemos ver que después de acercarnos repetidamente al punto de corte con el eje X mediante Zoom In y la mano Pan, la raíz buscada está comprendida entre 0.735 y 0.74.
Raíces simples Podemos utilizar el puntero del ratón para buscar los puntos en los que la función corta al eje X, añadiendo el comando ginput al final del script mouse_raiz x=0:0.1:pi/2; y=x-cos(x); plot(x,y) grid on xlabel('x') ylabel('y') title('x-cos(x)') [xRaiz,yRaiz] = ginput
En la ventana de comandos corremos el script mouse_raiz . La raíz buscada se señala mediante el cursor en forma de cruz de color negro. Se pulsa el botón izquierdo del ratón y luego la tecla Retorno. >> mouse_raiz xRaiz = 0.7379 yRaiz = -0.0022
Raíces múltiples Si la función tiene más de una raíz, se señala cada uno de los puntos de intersección de la función con el eje X con el puntero del ratón y se pulsa el botón izquierdo del ratón. Finalmente, se pulsa Retorno y aparecen las abscisas en el vector xRaiz y las ordenadas en el vector yRaiz de dichos puntos. Nos aproximaremos a las raíces que ya conocemos de una función para ilustrar este interesante procedimiento gráfico Guardamos el script mouse_raiz mediante File/Save As.. con el nombre mouse_raiz _1. Modificamos el script para representar la función coseno con el eje horizontal en grados y lo guardamos mediante File/Save o pulsando el icono correspondiente. Las dos raíces de cos( x)=0 en el intervalo (0,360) son 90° y 270°. x=0:1:360; y=cosd(x); plot(x,y)
grid on xlabel('x') ylabel('y') title('cos(x)') [xRaiz,yRaiz] = ginput
En la ventana de comandos corremos el script mouse_raiz _1. En la ventana gráfica, pulsamos dos veces el botón izquierdo del ratón: 1. cuando el puntero está situado aproximadamente en el punto de coordendas (90,0), pulsamos el botón izquierdo del ratón 2. cuando el puntero está situado en el punto (270,0), pulsamos el botón izquierdo del ratón 3. pulsamos Retorno, y obtenemos lo siguiente: >> mouse_raiz_1 xRaiz = 89.8618 yRaiz = -0.0088
270.5069 -0.0029
El procedimiento gráfico nos da una estimación de las raíces (puntos de intersección en la representación gráfica de la función con el eje X) de una ecuación trascendente en un intervalo ( a,b) dado.
La función MATLAB fzero La función fzero puede encontrar la raíz de una ecuación trascendente f ( x)=0. Su sintaxis es fzero( funcion, x0) Donde funcion es el nombre de la función cuyas raíces queremos determinar y x0 es el intervalo [a b] donde la función cambia de signo, es decir, el signo de f (a) es distinto al signo de f (b). x0 puede ser también un valor cercano a la raíz es decir, una primera aproximación. Podemos definir una función anónima y guardarla en el manejador func.Le pasamos la función anónima func a fzero. En la ventana de comandos, definimos la función y buscamos un intervalo [0.4 1] donde la función cambia de signo. Alternativamente, probamos una primara aproximación a la raíz buscada x0=0.6 >> func=@(x) cos(x)-x; >> func(0.4) ans = 0.5211 >> func(1) ans = -0.4597 >> r=fzero(func,[0.4 1]) %intervalo donde se encuentra la raíz r= 0.7391 >> fzero(func,0.6) %aproximación inicial a la raíz ans = 0.7391
Podemos definir explícitamente la función func y guardarla en el fichero func.m function y=func(x) y=cos(x)-x; end
Llamamos a fzero anteponiendo al nombre de la función el símbolo @ (véase al final de la página Funciones) >> fzero(@func,0.6) ans= 0.7391
A veces la función func definida y guardada en un fichero .M precisa de los valores de ciertos parámetros a, b function y=func(x,a,b) %código de la función end
A fzero solamente le podemos pasar el manejador (handle) de la función de variable x. La forma en que la función func conoce el valor de sus parámetros a y b es crear una función anónima f 1 en términos de func y se le pasamos su manejador a la función fzero f1=@(x)func(x,a,b);
p=fzero(f1,x0);
En el ejercicio "La ecuación de van der Waals" utilizaremos esta opción.
Funciones implícitas La siguiente función se estudia en la asignatura Energía Solar Fotovoltaica (tercer curso).
Donde x representa la diferencia de potencial e y la intensidad de la corriente eléctrica. Representamos la función f ( x,y) el intervalo [0,0.7] en el eje X y en el intervalo [0,3.1] en el eje Y. Utilizamos la función MATLAB contour que habitualmente se emplea para dibujar curvas de nivel. xx=linspace(0,0.7,50); yy=linspace(0,3.1,100); [x,y] = meshgrid(xx,yy); f=3-5e-12*(exp((x+y*0.01)/0.025)-1)-(x+y*0.01)/1000-y; contour(x,y,f,[0,0],'r'); xlabel('d.d.p. V'); ylabel('Intensidad'); title('Fotovoltaica') grid on
Una representación gráfica similar la obtenemos con la función ezplot Vamos a estudiar con más detalle esta función. En primer lugar, sabemos que para el valor de x próximo a 0.6 V (diferencia de potencial) la intensidad y vale cero. Empleamos la función MATLAB fzero para obtener un valor más preciso. f=@(x,y) 3-5e-12*(exp((x+y*0.01)/0.025)-1)-(x+y*0.01)/1000-y; g=@(x) f(x,0); Vmax=fzero(g,0.6); fprintf('Valor máximo de la d.d.p. %1.4f\n',Vmax)
En la ventana de comandos vemos este valor Valor máximo de la d.d.p. 0.6780
Dado un valor de la diferencia de potencial x en el intervalo [0,0.678], queremos obtener el correspondiente valor de la intensidad de la corriente y. Suministramos a fzero un valor próximo a la raíz buscada, por ejemplo, 3. Añadimos las siguientes líneas de código al script. .... x=0.6; g=@(y) f(x,y); y=fzero(g,3); fprintf('Para una d.d.p. de %1.4f le corresponde una intensidad de %1.4f\n',x, y)
En la ventana de comandos vemos este valor Para una d.d.p. de 0.6000 le corresponde una intensidad de 2.6214
Dado un vector de valores de la diferencia de potencial x el intervalo [0,0.678], queremos obtener el correspondiente vector de intensidades y. Finalmente, con estos dos vectores ( x,y) dibujaremos la curva diferencia de potencial - intensidad mediante el comando plot . El script completo es el siguiente f=@(x,y) 3-5e-12*(exp((x+y*0.01)/0.025)-1)-(x+y*0.01)/1000-y; g=@(x) f(x,0); Vmax=fzero(g,0.6); %intervalo de valores de la d.d.p. [0,Vmax] x=linspace(0,Vmax,100); n=length(x); opt=optimset('display','off'); %evita que fzero emita mensajes de aviso y=zeros(size(x)); for i=1:n y(i)=fzero(@(y) f(x(i),y),3,opt); end plot(x,y,'r') ylim([0 3.1]) xlabel('d.d.p. V') ylabel('Intensidad') title('Fotovoltaica')
Obtenemos una representación gráfica similar a la obtenida empleando la función MATLAB contour .
Sistemas de ecuaciones no lineales Calcular las raíces del sistema de dos ecuaciones: 2 f 1( x,y)=2 x - xy-5 x-1=0 2 f 2( x,y)= x+3log10 x- y =0
En la figura vemos el punto de intersección y mediante Tools/Data Cursor sus coordendas aproximadas.. x=linspace(2.5,5.5,50); y1=(2*x.^2-5*x-1)./x; y2=sqrt(x+3*log10(x)); plot(x,y1,'b',x,y2,'r')
Vamos a utilizar la función fsolve de MATLAB para obtener el punto de intersección. En primer lugar, tenemos que definir la función denominada sis_ecuaciones y guardarla en el correspondiente fichero . M. A esta función le tenemos que pasar los valores de x e y en el vector xn, y nos devuelve los valores de las funciones f 1( x,y) y f 2( x,y) en el vector s. function s=sis_ecuaciones(xn) x=xn(1); y=xn(2); s=[2*x^2-x*y-5*x-1, x*3*log10(x)-y^2]; end
Creamos el script sis_ecuaciones _ script para llamar al procedimiento numérico implementado en la función fsolve de MATLAB. x0 =[3 2]; %valor inicial [x,fval] = fsolve(@sis_ecuaciones,x0); fprintf('La solución es x=%1.3f, y=%1.3f\n',x(1),x(2)); fprintf('Valores de la función = %g\n',fval)
A fsolve tenemos que pasarle la función que hemos definido, sis_ecuaciones y la aproximaxión inicial ( x0, y0) que hemos encontrado anteriormente de forma gráfica. Esta función nos devuelve las coordenadas ( x,y) del punto de intersección buscado y los valores de las funciones f 1( x,y) y f 2( x,y) en dicho punto. En la ventana de comandos corremos el script sis_ecuaciones _ script >> sis_ecuaciones_script La solución es x=3.958, y=2.664 Valores de la función = -1.33692e-010 Valores de la función = -3.07305e-009
Energías Renovables ©EUITI de Eibar
