Método de Euler

La idea del método de Euler es muy sencilla y está basada en el significado geométrico de la derivada de una función en un punto dado.
Supongamos que tuviéramos la curva solución de la ecuación diferencial y trazamos la recta tangente a la curva en el punto dado por la condición inicial.

Debido a que la recta tangente aproxima a la curva en valores cercanos al punto de tangencia, podemos tomar el valor de la recta tangente en el punto

como una aproximación al valor deseado

Así, calculemos la ecuación de la recta tangente a la curva solución de la ecuación diferencial dada en el punto
De los cursos de Geometría Analítica, sabemos que la ecuación de la recta es:

donde m es la pendiente. En este caso, sabemos que la pendiente de la recta tangente se calcula con la derivada:

Por lo tanto, la ecuación de la recta tangente es :

Ahora bien, suponemos que es un punto cercano a y por lo tanto estará dado como
De esta forma, tenemos la siguiente aproximación:

De aquí, tenemos nuestra fórmula de aproximación:

Esta aproximación puede ser suficientemente buena, si el valor de h es realmente pequeño, digamos de una décima ó menos. Pero si el valor de h es más grande, entonces podemos cometer mucho error al aplicar dicha fórmula. Una forma de reducir el error y obtener de hecho un método iterativo, es dividir la distancia
en n partes iguales(procurando que estas partes sean de longitud suficientemente pequeña) y obtener entonces la aproximación en n pasos, aplicando la fórmula anterior n veces de un paso a otro, con la nueva h igual a
En una gráfica, tenemos lo siguiente:

Ahora bien, sabemos que: Para obtener únicamente hay que pensar que ahora el papel de lo toma el punto y por lo tanto, si sustituímos los datos adecuadamente, obtendremos que:

De aquí se ve claramente que la fórmula recursiva general, está dada por:

Esta es la conocida fórmula de Euler que se usa para aproximar el valor de

aplicándola sucesivamente desde
hasta

en pasos de longitud h.

Programa en matlab del método de Euler

function f
fprintf('\n \tRESOLUCION DE ECUACIONES DIFERENCIALES POR MEDIO METODO DE EULER\n')
f=input('\nIngrese la ecuacion diferencial de la forma: dy/dx=f(x,y)\n','s');
x0=input('\nIngrese el primer punto x0:\n');
x1=input('\nIngrese el segundo punto x1:\n');
y0=input('\nIngrese la condicion inicial y(x0):\n');
n=input('\nIngrese el numero de pasos n:\n');
h=(x1-x0)/n;
xs=x0:h:x1;
y1=y0;
fprintf('\n''it x0 x1 y1');
for i=1:n
it=i-1;
x0=xs(i);
x=x0;
x1=xs(i+1);
y=y0;
y1=y0+h*eval(f);
fprintf('\n%2.0f%10.6f%10.6f%10.6f\n',it,x0,x1,y1);
y0=y1;
end
fprintf('\n El punto aproximado y(x1) es = %10.6f\n',y1);

Solución
Respuesta

>> euler

RESOLUCION DE ECUACIONES DIFERENCIALES POR MEDIO METODO DE EULER

Ingrese la ecuacion diferencial de la forma: dy/dx=f(x,y)
sqrt(x^2+y^2)

Ingrese el primer punto x0:
2

Ingrese el segundo punto x1:
2.3

Ingrese la condicion inicial y(x0):
0.5

Ingrese el numero de pasos n:
3

'it x0 x1 y1
0 2.000000 2.100000 0.706155

1 2.100000 2.200000 0.927710

2 2.200000 2.300000 1.166470

El punto aproximado y(x1) es = 1.166470
>>

Método de Runge-Kutta

Programa en MatLab de Runge-Kutta de orden dos.

function f
fprintf('\n \tRESOLUCION DE ECUACIONES DIFERENCIALES POR MEDIO RUNGE-KUTTA DE ORDEN 4\n')
f=input('\n Ingrese la ecuacion diferencial dy/dx=\n','s');
x0=input('\n Ingrese el primer punto x0:\n');
x1=input('\n Ingrese el segundo punto x1:\n');
y0=input('\n Ingrese la condicion inicial y(x0):\n');
n=input('\n Ingrese el numero de pasos n:\n');
h=(x1-x0)/n;
xs=x0:h:x1;
fprintf('\n''it x0 y(x1)');
for i=1:n
it=i-1;
x0=xs(i);
x=x0;
y=y0;
k1=h*eval(f);
x=xs(i+1);
y=y0+k1;
k2=h*eval(f);
y0=y0+(k1+k2)/2;
fprintf('\n%2.0f%10.6f%10.6f\n',it,x0,y0);
end
fprintf('\n El punto aproximado y(x1) es = %8.6f\n',y0);


Programa de Runge-Kutta de orden cuatro.

function f
fprintf('\n \tRESOLUCION DE ECUACIONES DIFERENCIALES POR MEDIO RUNGE-KUTTA DE ORDEN 4\n')
f=input('\n Ingrese la ecuacion diferencial\n','s');
x0=input('\n Ingrese el primer punto x0:\n');
x1=input('\n Ingrese el segundo punto x1:\n');
y0=input('\n Ingrese la condicion inicial y(x0):\n');
n=input('\n Ingrese el numero de pasos n:\n');
h=(x1-x0)/n;
xs=x0:h:x1;
fprintf('\n''it x0 y(x1)');
for i=1:n
it=i-1;
x0=xs(i);
x=x0;
y=y0;
k1=h*eval(f);
x=x0+h/2;
y=y0+k1/2;
k2=h*eval(f);
x=x0+h/2;
y=y0+k2/2;
k3=h*eval(f);
x=x0+h;
y=y0+k3;
k4=h*eval(f);
y0=y0+(k1+2*k2+2*k3+k4)/6;
fprintf('\n%2.0f%10.6f%10.6f\n',it,x0,y0);
end
fprintf('\n El punto aproximado y(x1) es = %8.6f\n',y0);

Solucion

2. Respuesta

>> runge2

RESOLUCION DE ECUACIONES DIFERENCIALES POR MEDIO RUNGE-KUTTA DE ORDEN 4

Ingrese la ecuacion diferencial dy/dx=
4*exp(0.8*x)-0.5*y

Ingrese el primer punto x0:
0

Ingrese el segundo punto x1:
4

Ingrese la condicion inicial y(x0):
2

Ingrese el numero de pasos n:
4

'it x0 y(x1)
0 0.000000 6.701082

1 1.000000 16.319782

2 2.000000 37.199249

3 3.000000 83.337767

El punto aproximado y(x1) es = 83.337767

respuesta
>> runge4

RESOLUCION DE ECUACIONES DIFERENCIALES POR MEDIO RUNGE-KUTTA DE ORDEN 4

Ingrese la ecuacion diferencial
-2*x^3+12*x^2-20*x+8.5

Ingrese el primer punto x0:
0

Ingrese el segundo punto x1:
0.5

Ingrese la condicion inicial y(x0):
1

Ingrese el número de pasos n:
5

'it x0 y(x1)
0 0.000000 1.753950

1 0.100000 2.331200

2 0.200000 2.753950

3 0.300000 3.043200

4 0.400000 3.218750

El punto aproximado y(x1) es = 3.218750

Metodo de la División Sintética

function divisionsintetica
fprintf('Metodo de la Division Sintetica:\n');
p=input('Ingrese los Coeficientes del Polinomio como Vector:\n');
x0=input('Ingrese el numero a analizar:\n');
n=length(p);
a=p(1);
for j=1:n-1;
a(j+1)=p(j+1)+x0*a(j);
end
if(a(n)==0)
fprintf('ES UNA RAIZ=%5.3f\n',x0);
else
fprintf('NO ES UNA RAIZ\n');
end

Ejemplo:

Metodo de la Division Sintetica:
Ingrese los Coeficientes del Polinomio como Vector:
[1,3,-4]
Ingrese el numero a analizar:
-1
NO ES UNA RAIZ

Metodo de la Division Sintetica:
Ingrese los Coeficientes del Polinomio como Vector:
[1,3,-4]
Ingrese el numero a analizar:
1
ES UNA RAIZ=1.000

Metodo de la Biseccion

fprintf('\t\tmetodo de biseccion\n');
fun=input('ingrese la funcion:','s');
a=input('ingrese el punto inicial:');
b=input('ingrese el punto final:');
tol=input('ingrese la tolerancia:');
it=0;
x=a;
f_a=eval(fun);
x=b;
f_b=eval(fun);
if(f_a*f_b>0)
fprintf('\t\tEn este intervalo no hay raices\n');
fprintf('\t\tBuscar otro intervalo\n');
return
else
n=ceil(log((b-a)/tol)/log(2));
fprintf('\n=%10.6f\n',n);
fprintf(' it a b c f(a) f(b) f(c) abs(b-a)/2\n');
while 1
it=it+1;
c=(a+b)/2;
x=c;
f_c=eval(fun);
fprintf('%5.0f %10.9f %10.9f %10.9f %10.9f %10.9f %10.9f %10.9f\n',it,a,b,c,f_a,f_b,f_c,abs(b-a)/2);
if(abs(b-a)/2<=tol)
fprintf('n se satisface la tolerancia\n');
break end if(it>n)
fprintf('\numero de iteraciones excedido');
break
end
if(f_c*f_b<=0) a=c;
f_a=f_c;
else b=c;f_b=f_c;
end
end
fprintf('La raiz pedida es= %10.6f\n',c)
fprintf('La tolerancia es= %10.6f\n',abs(b-a)/2);
x=-100:0.1:100;
ezplot(fun),
grid on end


Ejemplo:

metodo de biseccion
ingrese la funcion:
3*x-2+exp(x)-x^2
ingrese el punto inicial:
0
ingrese el punto final:
1
ingrese la tolerancia:
10^(-5)

= 17.000000

n se satisface la tolerancia

La raiz pedida es= 0.257530

La tolerancia es= 0.000008

Metodo del Punto Fijo

function puntofijo
global fun
fprintf('Método del punto fijo:\n');
fun=input('Ingrese la función:\n','s');
x0=input('Ingrese el punto inicial:\n');
n=input('Ingrese el numero de iteraciones:\n');
it=0;
fprintf(' it x0 x1 x0-x1');
while(it it=it+1;
x=x0;
x1=eval(fun);
fprintf('\n%3.0f%15.10f%15.10f%15.10f\n',it,x0,x1,abs(x1-x0));
x0=x1;
end
fprintf('\n el punto fijo aproximado es=%10.6f\n',x1);
clf
hold on
ezplot('x');
ezplot fun
legend('y=x','fun')
hold off


Ejemplo:

Método del punto fijo:
Ingrese la función:
x^4+2*x^2-x-3
Ingrese el punto inicial:
1
Ingrese el numero de iteraciones:
10



el punto fijo aproximado es= 1.000000

Metodo de Newton Rapson

function newton
global fun dfun
fprintf('metodo de newton:\n');
fun=input('ingrese la funcion:\n','s');
x0=input('ingrese el punto inicial:\n');
tol=input('ingrese la tol:\n');
dfun=diff(fun);
it=0;
fprintf(' it x0 x1 x0-x1');
while(it<50) it="it+1;" x="x0;" x1="x0-(eval(fun)/eval(dfun));">fprintf('el procedimiento se completo satisfactoriamente:\n');
break
end
x0=x1;
end
fprintf('la raiz buscada es=%15.9f\n',x1);
ezplot(fun),grid on

Ejemplo:

metodo de newton:

ingrese la funcion:

exp(-x)-log(x)

ingrese el punto inicial:

1

ingrese la tol:

0.01

el procedimiento se completo satisfactoriamente:
la raiz buscada es= 1.309799389

Metodo de la Secante

function secante
global fun
fprintf('Metodo de la secante:\n');
fun=input('Ingrese la funcion:\n','s');
x0=input('Ingrese el primer punto inicial:\n');
x1=input('Ingrese el segundo punto inicial:\n');
tol=input('Ingrese la tol:\n');
it=0;
fprintf('it x0 x1 x2 x1-x2');
while(it<50) it="it+1;" x="x0;" f0="eval(fun);" x="x1;" f1="eval(fun);" x2="(x0*f1-x1*f0)/(f1-f0);"> fprintf('el procedimiento se completo satisfactoriamente:\n');
break
end
x0=x1;
x1=x2;
end
fprintf('la raiz buscada es=%15.9f\n',x2);
ezplot(fun),
grid on

Ejemplo:



Metodo de la secante:

Ingrese la funcion:

x^2+10*cos(x)

Ingrese el primer punto inicial:

1.5

Ingrese el segundo punto inicial:

2

Ingrese la tol:

10^(-5)

el procedimiento se completo satisfactoriamente:

la raiz buscada es= 1.968872938

Metodo de Falsa Posicion

function falsaposicion
global fun
fprintf('\n metodo de la falsa posicion\n');
fun=input('\nIngrese la funcion:\n','s');
x0=input('\nIngrese el valor inferior del intervalo\n');
x1=input('\nIngrese el valor superior del intervalo\n');
tol=input('\nIngrese la tolerancia\n')
fprintf('\n it x0 x1 x2 f0 f1 f2 abs(f2)\n');
f2=1;
it=0;
x=x0;
f0=eval(fun);
x=x1;
f1=eval(fun);
if(f0*f1>0)
fprintf('\n En este intervalo no hay raices\n');
fprintf('\n Buscar otro intervalo \n');
return
end
while abs(f2)>tol
it=it+1;
x2=(x0*f1-x1*f0)/(f1-f0);
x=x2;
f2=eval(fun);
fprintf('\n%3.0f%12.7f%12.7f%12.7f%12.7f%12.7f%12.7f%12.7f\n',it,x0,x1,x2,f0,f1,f2,abs(f2));
if(f0*f2<0) x1="x2;f1=" x0="x2;f0=" es =" %10.7f\n',x2);" face="arial" color="#cc0000">


Ejemplo:

metodo de la falsa posicion

Ingrese la funcion:
3*x-2+exp(x)-x^2

Ingrese el valor inferior del intervalo
0
Ingrese el valor superior del intervalo
1
Ingrese la tolerancia
10^(-5)
tol =
1.0000e-005

it	x0	x1	x2	f0	f1	f2	abs(f2)
1	0.0000000	1.0000000	0.2689414	-1.0000000	2.7182818	0.0430733	0.0430733
2	0.0000000	0.2689414	0.2578356	-1.0000000	0.0430733	0.0011537	0.0011537
3	0.0000000	0.2578356	0.2575385	-1.0000000	0.0011537	0.0000310	0.0000310
4	0.0000000	0.2575385	0.2575305	-1.0000000	0.0000310	0.0000008	0.0000008

La raiz buscada es = 0.2575305

Metodo de Newton Modificado

function newtonmodificado
global fun dfun ddfun
fprintf('metodo de newtonmodificado:\n');
fun=input('ingrese la funcion:\n','s');
x0=input('ingrese el punto inicial:\n');
tol=input('ingrese la tolerancia:\n');
dfun=diff(fun);
ddfun=diff(diff(fun));
it=0;
fprintf('it x0 x1 x0-x1');
while(it<50) it="it+1;" x="x0;" x1="x0-((eval(fun)*eval(dfun))/((eval((dfun)^2))-(eval(fun)*eval(ddfun))));"> fprintf('el procedimiento se completo satisfactoriamente:\n');
break
end
x0=x1;
end
fprintf('la raiz buscada es =%15.9f\n',x1);
ezplot(fun),grid on


Ejemplo:


metodo de newtonmodificado:
ingrese la funcion:
exp(x)-x-1
ingrese el punto inicial:
0.5
ingrese la tolerancia:
10^(-5)

it	x0	x1	x0-x1
1	0.500000000	-0.049299708	0.549299708
2	-0.049299708	-0.000398480	0.048901228
3	-0.000398480	-0.000000026	0.000398453
4	-0.000000026	-0.000000014	0.000000012

el procedimiento se completo satisfactoriamente:

la raiz buscada es = -0.000000014

Metodo de los Minimos Cuadrados

function correl
xi=input('ingrese los valores de xi:\n');
yi=input('ingreso los valores de yi:\n');
xm=input('ingreso un valor de x para analisar la pendiente:\n');
n=length(xi);
p=xi.';
q=yi.';
r=p.*q;
s=p.*p;
a1=sum(xi);
a2=sum(yi);
a3=sum(r);
a4=sum(s);
m=a2/n;
j=sum((yi-m).^2);
A=[a4 a1;a1 n];
B=[a3 a2];
X=inv(A)*(B');
l=X(1);r=X(2);
fprintf('\nLos coeficientes de la ecuacion son:\n')
fprintf(' a b\n')
fprintf('%8.8f%10.8f\n',l,r)
fun=input('\nLa funcion lineal es y=','s');
y=l*xi+r;
t=sum((yi-y).^2);
r2=1-(t/j);
z=(r2)^0.5;
m=diff(fun);
x=xm;
n=eval(m);
fprintf('\n El coeficiente de determinacion r^2=');
fprintf('%5.5f\n',r2);
if(n>0)
fprintf('\n El coeficiente de correlacion r=');
fprintf('%5.5f\n',z);
elsefprintf('\n El coeficiente de correlacion r=');
fprintf('%5.5f\n',-z);
end
hold on
ezplot(fun),grid on,
plot(xi,yi,'m.'),
hold off


Ejemplo: Funcion Lineal

ingrese los valores de xi:

[4 5 2 5 6 7 1 8 3 7]

ingreso los valores de yi:

[5 6 4 5 7 10 3 11 4 9]

ingreso un valor de x para analizar la pendiente:

2
Los coeficientes de la ecuacion son:

a b

1.13025210 0.97478992

La funcion lineal es y=1.13025210*x+0.97478992

El coeficiente de determinacion r^2=0.88900

El coeficiente de correlacion r=0.94287

function correlog2
xi=input('ingrese los valores de xi:\n');
yi=input('ingreso los valores de yi:\n');
n=length(xi);
p=log(-xi);
q=yi;
r=p.*q;
s=p.*p;
a1=sum(p);
a2=sum(q);
a3=sum(r);
a4=sum(s);
m=sum(yi)/n;
j=sum((yi-m).^2);
A=[a4 a1;a1 n];
B=[a3 a2];
X=inv(A)*(B');
l=X(1);r=X(2);
fprintf('\nLos coeficientes de la funcion son:\n')
fprintf(' a b\n')
fprintf('%8.12f%10.12f\n',l,r)
y=input('\nLa funcion logaritmica natural es y=','s');
h=l*log(-xi)+r;
z=sum((yi-h).^2);
r2=1-(z/j);
k=(r2)^0.5;
fprintf('\n El coeficiente de determinacion r^2=');
fprintf('%5.5f\n',r2);
fprintf('\n El coeficiente de correlacion r=');
fprintf('%5.5f\n',k);
hold on
ezplot(y),grid on,
plot(xi,yi,'m.'),
hold off