1-MÉTODO DE LA BISECCIÓN

El método de la bisección nos permitirá resolver funciones no analíticas mediante aproximaciones. Para ello habrá que escoger un intervalo en el cual la función deberá cambiar de signo, cumpliendo lo siguiente:

Intervalo [a,b] ! f(a)*f(b)<=0

a b

Cogiendo el valor medio © de [a,b], tengo que comprobar:

f(a)*f(c)<=0 ! entonces realizaremos la siguiente asignación b!c

sino a!c

Repetiremos recursivamente hasta que el resultado tenga un error menor al margen que nosotros deseamos, siendo la raíz el valor medio.

Si quisiéramos calcular la bisección de la función:

F(x)=-19 (X-0.5) (X-1) + e^(x) - e^(-2x)

En primer lugar realizaremos un rastreo de la función para observar los cambios de signo, con lo cual sabremos en que intervalos debemos realizar la bisección.

Realizaremos el rastreo mediante el siguiente programa:

program intervalo;

var

i:integer;

f:real;

begin

for i:=-10 to 10 do

begin

f:=-19*(i-0.5)*(i-1)+exp(i)-exp(-2*i);

writeln ('X= ',i,' Y= ',f);

end;

end.

X	Y
-10	-4.8516738991E+08
...	...
0	-9.5000000000E+00
1	2.5829465452E+00
2	-2.1129259540E+01
...	...
6	-1.1907121265E+02
7	3.5563315759E+02
...	...
10	2.0401965795E+04

Analizando los resultados obtenidos observamos que hay tres cambios de signo, por lo que deducimos que tendremos que realizar la bisección en tres zonas , las cuales son:

1ª bisección ! [0,1]

2ª bisección ! [1,2]

3ª bisección ! [6,7]

Debemos tener en cuenta para que se cumpla la condición (f(c)*f(a)<=0) que deberá haber un número impar de raíces dentro del intervalo definido, si no podría haber irregularidades tales como: Discontinuidades o zonas de ruptura, que nos causarían problemas de cálculo.

Una vez rastreada la zona, para realizar la bisección utilizaremos el siguiente programa:

program biseccion; const tolerancia=0.0000001; var a,b,c:real;

function f (x:real):real;

begin

f:=-19*(x-0.5)*(x-1)+exp(x)-exp(-2* x);

end;

begin

writeln ('Introduce los valores del intervalo:');

readln (a,b);

repeat

c:=(a+b)/2;

if f(c)*f(a)<=0 then

b:=c

else

a:=c;

until abs(b-a)<tolerancia;

writeln ('La raiz es: ',(b+a)/2);

end.

Los resultados obtenidos en cada bisección son:

[0,1] ! RAIZ= 4.0626206994 E-01

[1,2] ! RAIZ= 1.2398312986 E+00

[6,7] ! RAIZ= 6.4089687765 E+00

Para definir una variedad de funciones, utilizaremos como ejemplo la siguiente función genérica, en la cual variaran los valores de las constantes:

f(x)= C1*X^(C2)+C3*X^(C4)+C5*e^(C6*X)

Para calcular las raíces de la función genérica, utilizaremos un programa en el cual, una vez introducidas las constantes realizará el rastreo en el intervalo definido. Una vez analizados estos resultados, aplicará la bisección y nos dará la raiz.

program generico; var c1,c2,c3,c4,c5,c6,a,b:integer; function pot (a:integer;b:real):real; var i:integer;

potencia:real;

begin

potencia:=1;

for i:=1 to a do

potencia:=potencia*b;

pot:=potencia;

end;

procedure rastreo (c1,c2,c3,c4,c5,c6:integer);

var

i,a,b:integer;

f:real;

begin

writeln ('Introduce el intervalo para realizar el rastreo:');

readln (a,b);

for i:=a to b do

begin

f:=c1*pot(c2,i)+c3*pot(c4,i)+c5*exp(c6*i);

writeln ('X= ',i,' Y= ',f);

end;

end;

procedure raiz (c1,c2,c3,c4,c5,c6:integer);

var

c,a,b:real;

begin

writeln ('Introduce el intervalo a observar:');

readln (a,b);

repeat

c:=(a+b)/2;

if (c1*pot(c2,c)+c3*pot(c4,c)+c5*exp(c6*c))*(c1*pot(c2,a)+

c3*pot(c4,a)+c5*exp(c6*a))<=0 then b:=c else a:=c;

until abs(b-a)<0.0000001; writeln ('La solucion es: ',(b+a)/2); end; begin writeln ('Introduce las constantes:');

readln (c1,c2,c3,c4,c5,c6);

rastreo (c1,c2,c3,c4,c5,c6);

raiz(c1,c2,c3,c4,c5,c6);

end.

Si utilizamos un intervalo de rastreo de [-10,10] y definimos diferentes constantes, con los valores obtenidos rellenaremos la gráfica contigua:

C1	C2	C3	C4	C5	C6	INTERVALO	RAIZ
2	5	3	0	1	0	[-2,-1]	-1.1486983597 E+00
-1	2	3	-4	5	-1	[1,2]	1.9301833808 E+00
1	1	1	0	1	0	[-3,-2]	-2.0000000298 E+00

Con la función genérica definida, introduciéndoles valores diferentes a las constantes, podremos definir diferentes funciones, siendo siempre de la familia de las potencias y exponenciales.

2-MÉTODO DEL PUNTO FIJO

Partiendo de una función f(x)=0, despejando x obtendremos otra función x=g(x). La solución de nuestro problema será donde se corten dicha función y la función y=x.

Para llegar a la solución deberemos introducir una condición de partida (Xo).

g(x) y=x

 ! solución

Para obtener el resultado mediante el método del punto fijo, comprobaremos que g(x) es derivable en el intervalo [a,b].

Una vez realizada la derivada, para comprobar si la función converge o diverge, realizaremos la siguiente sustitución: g'(Xo)

Si |g'(Xo)|<1 diremos que converge hacia , diferenciando dos casos:

y=x y=x

g(x)

g(x)

-1<g'(x)<=0 0<=g'(x)<1

Sin embargo,si |g'(x)|>1 diremos que diverge de , diferenciando los siguientes dos casos:

y=x y=x

g(x) g(x)

g'(x)<-1 g'(x)>1

Las condiciones anteriormente citadas, son solo suficientes y no necesarias, ya que puede suceder que |g'(x)|>1 y sea convergente.

A continuación resolveremos el siguiente ejercicio:

F(x)=x-e^(1/x)

Xo=1.5

Tolerancia=0.00005

A la hora de despejar X, hemos analizado las siguientes opciones:

1/ X=g(x)=e^(1/x) ! g´(x)=-e^(1/x)/(x*x) ! |g'(x)|<1! CONVERGENTE

2/ X=g(x)=1/ln(x) ! g'(x)=x ! |g´(x)|>1 ! DIVERGENTE

3/ X=g(x)=x/2+e^(1/x)/2 ! g´(x)=1/2-e^(1/x)/2*(x*x) ! |g´(x)|<1 ! CONVERGENTE

program punto_fijo; var x,tol,aux,x0:real; cont:integer; function calculo (x0:real):real; begin calculo:=exp(1/x);* end; begin writeln ('Introduce el valor inicial de X:');

readln (x0);

writeln ('Introduce la tolerancia deseada:');

readln (tol);

cont:=0;

repeat

x:=calculo (x0);

cont:=cont+1;

aux:=x0;

x0:=x;

until (abs(aux-x0)<tol);

writeln ('El calculo se ha realizado ',cont,' veces');

writeln ('La soluci¢n es = ',x0);

end.

• Se sustituirá por la función g(x)=x/2+(e^(1/x)/2), en el tercer caso

No hemos tenido en cuenta el segundo caso porque es divergente.

Las soluciones obtenidas mediante este programa son:

FUNCIONES	Nº DE INTERACCIONES	SOLUCIÓN
g(x)=e^(1/x)	18	1.7632117954
g(x)=1/ln x	DIVERGENTE
g(x)=x/2+(e^(1/x)/2)	7	1.7632189846

A la conclusión que hemos llegado al observar los resultados, es que la manera de despejar X tiene su importancia a la hora de calcular la solución. Basándonos en el ejemplo anterior podemos ver que la tercera opción es la más eficaz, ya que con menos interacciones obtenemos prácticamente la misma solución.

Definiremos una función general, mediante el cual, variando las constantes se podrán realizar diferentes cálculos:

g(x)=(C1*X^(C2)+C3*X^(C4)+C5*e^(C6*X))/C7

g´(x)=(C1*C2*X^(C2-1)+C3*C4*X^(C4-1)+(C5*e^(C6*X))*(C6+1))/C7

Para este cálculo utilizaremos el siguiente programa:

program punto_fijo;

var

x,tol,aux,x0:real;

cont,c1,c2,c3,c4,c5,c6,c7:integer;

function pot (a:real;b:integer):real;

var

potencia:real;

i:integer;

begin

potencia:=1;

for i:=1 to b do

potencia:=potencia*a;

pot:=potencia;

end;

function deribada (x0:real;c1,c2,c3,c4,c5,c6,c7:integer):boolean;

var

valor:real;

begin

valor:=(c1*c2*pot(x0,c2-1)+c3*c4*pot(x0,c4-1)+c5*exp(c6*x)*(c6+1))/c7;

if (abs(valor)<1) then

deribada:=true

else

deribada:=false;

end;

function calculo (x0:real;c1,c2,c3,c4,c5,c6,c7:integer):real;

begin

calculo:=(c1*pot(x0,c2)+c3*pot(x0,c4)+c5*exp(x0*c6))/c7;

end;

begin

writeln ('Introduce el valor inicial de X:');

readln (x0);

writeln ('Introduce la tolerancia deseada:');

readln (tol);

writeln ('Introduce las constantes:');

readln (c1,c2,c3,c4,c5,c6,c7);

if deribada (x0,c1,c2,c3,c4,c5,c6,c7) then

begin

cont:=0;

repeat

x:=calculo (x0,c1,c2,c3,c4,c5,c6,c7);

cont:=cont+1;

aux:=x0;

x0:=x;

until (abs(aux-x0)<tol);

writeln ('El calculo se ha realizado ',cont,' veces');

writeln ('La soluci¢n es = ',x0);

end

else

writeln ('La funci¢n es divergente');

end.

Los resultados obtenidos con las diferentes constantes:

C1	C2	C3	C4	C5	C6	C7	Nº INT	RESULTADO
1	2	-2	0	1	1	3	5	-0.39027166481
1	1	3	-1	2	0	4	12	1.6666665673
0	0	1	0	1	1	1	DIVERGENTE

3-MÉTODO DE NEWTON-RAPHSON

El método de Newton se obtiene a partir del desarrollo de Taylor:

Taylor ! f(X)=f(Xo)+f´(Xo)*(X-Xo)+(f´´(Xo)/2!)*(X-Xo)^2+...

f(X)=f(Xo)+f´(Xo)*(X-Xo)+0*h^2 ! h=X-Xo

Debemos seleccionar una cantidad de sumas dentro de una serie infinita, por lo tanto realizaremos un error de truncamiento.

Con este método se toma un valor de partida Xo y se van calculando valores sucesivos hasta que Xn este dentro de la tolerancia requerida, la cual consideraremos que es el resultado correcto.

Para ello la ecuación truncada de Taylor la igualamos con 0, llamandolo g(X). De esta despejaremos X obteniendo la ecuación que utilizaremos para el cálculo de la solución.

f(X)=0

g(X)=f´(Xo)(X-Xo)+f(Xo)

! !

m b

g(x)=m*x+b=0 ! Xn=X(n-1)-f(X(n-1))/f´(X(n-1))

La ventaja de este método con respecto a las anteriores es que con pocas interacciones obtendremos el resultado. Sin embargo, tiene la pega de necesitar la derivada de la función, para dicho cálculo, aveces nos será necesario emplear algun otro método númerico.

A continuación realizaremos varios cálculos utilizando el método de Newton-Raphson:

f(X)=X^3-155 ! f`(X)=3 X^2

Tolerancia=0.0000005

Realizaremos el cálculo para dos valores de Xo:

Xo=5 Xo=10

program newton;

var

x0,tol,x,aux:real;

cont:integer;

function calculo (x0:real):real;

var

result:real;

begin

result:=((x0*x0*x0)-155)/(3*(x0*x0));

calculo:=x0-result;

end;

begin

writeln ('Introduce el valor inicial:');

readln (x0);

writeln ('Introduce el valor de la tolerancia:');

readln (tol);

cont:=0;

repeat

x:=calculo(x0);

cont:=cont+1;

aux:=x0;

x0:=x;

until (abs(aux-x0)<tol);

writeln ('El numero de operaciones realizadas = ',cont);

writeln ('La soluci¢n es = ',x0);

end.

Los resultados obtenidos son los siguientes:

VALOR INICIAL (Xo)	Nº INTERACCIONES	RESULTADO
5	4	5.3716853549
10	6	5.3716853549

Observando el resultado obtenido, nos damos cuenta que no hay cambios para distintos valores iniciales.

f(X)=e^X-X-2 ! f´(X)=e^X-1

Tolerancia=0.0000005

Realizaremos el cálculo para dos valores de Xo:

Xo=-10 Xo=10

program newton;

var

x0,tol,x,aux:real;

cont:integer;

function calculo (x0:real):real;

var

result:real;

begin

result:=(exp(x0)-x0-2)/(exp(x0)-1);

calculo:=x0-result;

end;

begin

writeln ('Introduce el valor inicial:');

readln (x0);

writeln ('Introduce el valor de la tolerancia:');

readln (tol);

cont:=0;

repeat

x:=calculo(x0);

cont:=cont+1;

aux:=x0;

x0:=x;

until (abs(aux-x0)<tol);

writeln ('El numero de operaciones realizadas = ',cont);

writeln ('La soluci¢n es = ',x0);

end.

VALOR INICIAL (Xo)	Nº INTERACCIONES	RESULTADO
-10	4	-1.84056604
10	14	1.1461932206

En este caso vemos que con distintos valores iniciales cambia el resultado, esto es debido a la función, porque tiene dos raíces en el intervalo

[-10,10].

Para distintas alternativas, introduciremos una función genérica:

f(X)=C1*X^C2+C3*X^C4+C5*e^(X*C6)

f´(X)=C1*C2*X^(C2-1)+C3*C4*X^(C4-1)+C5*C6*e^(C6*X)

Tolerancia=0.000001

Valor inicial ! Xo=0

program newton;

var

x0,tol,x,aux:real;

cont,c1,c2,c3,c4,c5,c6:integer;

function pot (a:real;b:integer):real;

var

multi:real;

i:integer;

begin

multi:=1;

for i:=1 to b do

multi:=multi*a;

pot:=multi;

end;

function calculo (x0:real):real;

var

result:real;

begin

result:=(c1*pot(x0,c2)+c3*pot(x0,c4)+c5*exp(c6*x0))/

(c1*c2*pot(x0,c2- 1)+c3*c4*pot(x0,c4-1)+c5*c6*exp(c6*x));

calculo:=x0-result;

end;

begin

writeln ('Introduce el valor inicial:');

readln (x0);

writeln ('Introduce el valor de la tolerancia:');

readln (tol);

writeln ('Introduce el valor de las constantes:');

readln (c1,c2,c3,c4,c5,c6);

cont:=0;

repeat

x:=calculo(x0);

cont:=cont+1;

aux:=x0;

x0:=x;

until (abs(aux-x0)<tol);

writeln ('El n£mero de operaciones realizadas = ',cont);

writeln ('La soluci¢n es = ',x0);

end.

Para distintos valores que hemos introducido, los resultados han sido:

C1	C2	C3	C4	C5	C6	Nº INT	RESULTADO
1	2	-2	0	1	1	6	0.53727444917
1	2	-3	0	2	1	5	0.36104234240

Estos cálculos se podrían realizar introduciendo diferentes valores de las constantes. Teniendo cuidado que no nos quede el valor 0 en el divisor.

0

12