lab8--

Vista previa del material en texto

UNIVERSIDAD DE CONCEPCION
FACULTAD DE CIENCIAS
FISICAS Y MATEMATICAS
DEPARTAMENTO DE INGENIERIA MATEMATICA
Cálculo Numérico (521230)
Laboratorio 8
Ecuaciones no lineales
El objetivo de este laboratorio es aprender técnicas para la resolución numérica de ecuaciones y
sistemas de ecuaciones no lineales.
1. (a) Haga un programa que calcule la ráız de una ecuación f(x) = 0 mediante el método de bisección.
Los datos del programa deben ser la función f , los extremos del intervalo [a, b] donde se busca
la ráız, y la tolerancia tol con la que se desea calcular ésta.
El programa debe tener una salida de error en el caso en que la función f no cambie de signo
en el intervalo inicial.
(b) Calcule con el programa anterior todas las ráıces de las siguientes ecuaciones con error menor
que 10−4:
x2 = 2, x3 − 3x + 1 = 0 y cosx = x.
2. El archivo newton.m (bájelo de la página web del curso o solićıtelo al ayudante) contiene el siguiente
programa para el cálculo de la ráız de una ecuación f(x) = 0 mediante el método de Newton-
Raphson:
function raiz=newton(f,Df,x0,tol,maxit)
k=0;
raiz=x0;
corr=tol+1;
while (k<maxit) & (abs(corr)>tol)
k=k+1;
xk=raiz;
fxk=feval(f,xk);
Dfxk=feval(Df,xk);
if (Dfxk==0)
error(’La derivada de la funcion se anula.’)
end
corr=fxk/Dfxk;
raiz=xk-corr;
end
if (abs(corr)>tol)
error(’Se excedio el numero maximo de iteraciones.’)
end
1
(a) Calcule mediante este programa las ráıces de las ecuaciones del Ej. 1(b) con error menor que
10−12.
(b) Modifique el programa para que permita resolver sistemas de n ecuaciones no lineales con n
incógnitas.
Sugerencia: Dado que Matlab trabaja del mismo modo con variables numéricas o vecto-
riales, sólo hace falta modificar ligeramente las siguientes ĺıneas del programa para que éste
sirva para sistemas de ecuaciones:
ĺıneas 6 y 16: en lugar de abs(corr) debe usarse alguna medida del vector de correcciones;
ĺınea 9: para sistemas, el programa puede fallar aunque la diferencial Df no sea nula;
ĺınea 10: el mensaje de error debe adaptarse correspondientemente;
ĺınea 11: Df(x(k))−1f(x(k)) no se calcula mediante fxk/Dfxk en el caso vectorial;
(c) Calcule con el programa anterior todas las ráıces de los siguientes sistemas de ecuaciones con
error menor que 10−12:
{
x2 + xy + y2 = 1
y = x2
e
{
y = e−x
x = ey
3. El comando Matlab fzero permite calcular la ráız de una ecuación f(x) = 0 cercana a un punto
dado x0. Este comando combina de manera automática un método inicial de convergencia garan-
tizada con otro final de convergencia veloz.
(a) Utilice el help de Matlab para conocer la sintaxis del comando fzero.
(b) Calcule, mediante este comando, las ráıces de las ecuaciones del Ej. 1(b), primero con la
tolerancia prefijada por el comando y luego con error menor que 10−12.
4. Resuelva los siguientes problemas:
(a) Calcular el área encerrada por las siguientes curvas:
y = ex−x
2
e y = arctanx2.
(b) Calcular todos los valores de x para los que
∫ x
0
sen t2 dt =
√
π
4
.
RAD/RRS/RRA/MSC
http://www.ing-mat.udec.cl/pregrado/asignaturas/521230/
17/06/06
2
Soluciones propuestas
1. Ej. 1 (a). File bisec.m :
function raiz=bisec(funct,a,b,tol)
fa=feval(funct,a);
fb=feval(funct,b);
if (fb*fa>0)
error(’La funcion tiene el mismo signo en ambos extremos.’)
end
while (abs(b-a)>tol)
if (fa*fb==0)
if (fa==0)
raiz=a;
else
raiz=b;
end
else
raiz=(a+b)/2;
fraiz=feval(funct,raiz);
if (fa*fraiz>0)
a=raiz;
fa=fraiz;
else
b=raiz;
fb=fraiz;
end
end
end
3
2. Ej. 1 (b). Primero debe graficarse cada función (o eventualmente un par de funciones, como en
la tercera ecuación) para localizar las ráıces buscadas. Las funciones pueden definirse mediante el
comando inline, por ejemplo,
>> f2=inline(’x.^3-3*x+1’);
>> x=(-3:.1:3);
>> plot(x,feval(f2,x),x,zeros(length(x),1);
>> raiz1=bisec(f2,-2,-1,1.e-4)
raiz1 =
-1.8793
>> raiz2=bisec(f2,0,1,1.e-4)
raiz2 =
0.3474
>> raiz3=bisec(f2,1,2,1.e-4)
raiz3 =
1.5320
o bien mediante un programa function, por ejemplo,
File f3.m :
function y=f3(x)
y=cos(x)-x;
Salida:
>> raiz=bisec(’f3’,0,1,1.e-4)
raiz =
0.7391
Nótese que si la función se define mediante un programa function, su nombre debe pasarse entre
comillas como parámetro a otro programa o comando (por ejemplo, bisec(’f3’,0,1,1.e-6)).
En cambio, si se define inline, debe pasarse sin comillas (por ejemplo, bisec(f2,-2,-1,1.e-6)).
3. Ej. 2 (a). Por ejemplo:
>> f3=inline(’cos(x)-x’);
>> Df3=inline(’-sin(x)-1’);
>> tol=1.e-12;
>> maxit=10;
>> x0=1;
>> format long
>> raiz=newton(f3,Df3,x0,tol,maxit)
raiz =
0.73908513321516
4
4. Ej. 2 (b). File newton.m :
function raiz=newton(f,Df,x0,tol,maxit)
k=0;
raiz=x0;
corr=tol+1;
while (k<=maxit) & (norm(corr,inf)>tol)
k=k+1;
xk=raiz;
fxk=feval(f,xk);
Dfxk=feval(Df,xk);
if (rank(Dfxk)<length(Dfxk))
error(’La diferencial de la funcion es singular.’)
end
corr=Dfxk\fxk;
raiz=xk-corr;
end
if (norm(corr,inf)>tol)
error(’Se excedio el numero maximo de iteraciones.’)
end
5. Ej. 2 (c). Primero deben localizarse las raices para lo que debe esbozarse el gráfico de las respectivas
curvas. En este caso puede ser más fácil hacerlo anaĺıticamente que mediante Matlab. Al menos
será necesario realizar algún trabajo anaĺıtico antes de usar el comando plot de Matlab.
File Ej8 2c.m :
f1=inline(’[x(1)^2+x(1)*x(2)+x(2)^2-1; x(2)-x(1)^2]’);
Df1=inline(’[2*x(1)+x(2) x(1)+2*x(2); -2*x(1) 1]’);
tol=1.e-12;
maxit=10;
x0=[1;1];
raiz1_1=newton(f1,Df1,x0,tol,maxit)
x0=[-1;1];
raiz1_2=newton(f1,Df1,x0,tol,maxit)
f2=inline(’[x(2)-exp(-x(1)); x(1)-exp(x(2))]’);
Df2=inline(’[exp(-x(1)) 1; 1 -exp(x(2))]’);
x0=[1;1];
raiz2=newton(f2,Df2,x0,tol,maxit)
Salida:
>> Ej8_2c
raiz1_1 =
0.68232780382802
0.46557123187677
raiz1_2 =
-1
1
raiz2 =
1.30979958580415
0.26987413757345
5
6. Ej. 3 (b). Por ejemplo:
>> f1=inline(’x.^2-2’);
>> raiz1=fzero(f1,1.5)
raiz1 =
1.4142
>> format long
>> raiz2=fzero(f1,1.5,1.e-12)
raiz2 =
1.41421356237309
7. Ej. 4(a). Primero hay que dibujar las gráficas de las funciones, luego determinar los puntos de
intersección y finalmente integrar la diferencia de las dos funciones:
File Ej8 2c.m :
x=-2:.01:3;
plot(x,exp(x-x.^2),x,atan(x.^2))
f=inline(’exp(x-x.^2)-atan(x.^2)’);
a=fzero(f,-.5)
b=fzero(f,1)
integ=quad(f,a,b)
Salida:
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3
0
0.5
1
1.5
Salida:
>> Ej8_4a
Zero found in the interval: [-0.34, -0.66].
a =
-0.6196
Zero found in the interval: [0.88686, 1.1131].
b =
1.1080
integ =
1.2372
6
8. Ej. 4(b). Primero debe graficarse la función
F (x) =
∫ x
0
sen t2 dt −
√
π
4
para localizar sus ráıces. Luego se determinan éstas. Nótese que para evaluar la función F (x),
tanto para graficarla como para calcular sus ráıces, es necesario evaluar la integral.
File F.m :
function y=F(x)
n=length(x);
f=inline(’sin(t.^2)’);
for i=1:n
y(i)=quad(f,0,x(i))-sqrt(pi)/4;
end
>> x=-5:.1:5;
>> plot(x,F(x),x,zeros(length(x),1))
−5 −4 −3 −2 −1 0 1 2 3 4 5
−1.5
−1
−0.5
0
0.5
>> x=1:.025:3;
>> plot(x,F(x),x,zeros(length(x),1))
1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5
7
Salida:
>> raiz1=fzero(’F’,1)
Zero found in the interval: [0.84, 1.16].
raiz1 =
1.1459
>> raiz2=fzero(’F’,2.4)
Zero found in the interval: [2.3321, 2.4679].
raiz2 =
2.4347
>> raiz3=fzero(’F’,2.6)
Zero found in the interval: [2.5265, 2.6].
raiz3 =
2.5779
8