INFORME FINAL

MÉTODOS NUMÉRICOS PARA RESOLVER ECUACIONES DIFERENCIALES

PROFESOR: DOCTOR ANGEL SANGIACOMO CARAZAS

ALUMNO: LARICO ALVIZ JOEL NELSON

17 DE AGOSTO DE 2022
INFORME SOBRE LOS ALGORITMOS DE MÉTODOS NUMÉRICOS PARA LA
SOLUCIÓN DE ECUACIONES DIFERENCIALES

REDACTADO POR

Larico Alviz Joel Nelson

ROGRAMAS Y/O HERRAMIENTAS UTILIZADAS

 Octave
 Calculadora Científica
 Libro de MNPSED

COMPETENCIA ESPECIFICA

 Programar y estudiar los diferentes algoritmos acerca de los diferentes

métodos para facilitar el cálculo y generar sus gráficas en una Ecuación
Diferencial, con ayuda de programas como Matlab u Octave.

OBJETIVOS GENERALES

Al terminar este capítulo el lector estará en condiciones de:

Resolver ecuaciones diferenciales ordinarias analíticamente, usando la función
incorporada dsolve de Matlab.
Resolver ecuaciones diferenciales ordinarias con valor inicial usando la función
incorporada ode45 de Matlab.
Resolver ecuaciones diferenciales ordinarias con valor inicial usando un programa
codificado en Matlab aplicando el método de Euler simple.
Resolver ecuaciones diferenciales ordinarias con valor inicial usando un programa
codificado en Matlab aplicando el método de Euler modificado.
Resolver sistemas de ecuaciones diferenciales ordinarias con valor inicial
usando programas codificados en Matlab aplicando los métodos de Euler
simple y modificado.
Resolver ecuaciones diferenciales ordinarias con valor inicial usando un programa
codificado en Matlab aplicando el método de Runge-Kutta de cuarto orden.
Resolver sistemas de ecuaciones diferenciales ordinarias con valor inicial
usando
programas codificados en Matlab aplicando los métodos de Runge-Kutta de cuarto.
Resolver ecuaciones diferenciales rígidas usando las funciones incorporadas ode15i
 ALGORITMOS APLICADOS EN CLASE

MÉTODO DE EULER

¿Qué es el Método de Euler?

En matemática y computación, el método de Euler, llamado así en honor a Leonhard

Euler, es un procedimiento de integración numérica para resolver ecuaciones
diferenciales ordinarias (EDO) a partir de un valor inicial dado. El método de Euler
es el más simple de los métodos numéricos para resolver un problema de valor
inicial, y el más simple de los Métodos de Runge-Kutta. El método de Euler es
nombrado por Leonhard Euler, quien lo trató en su libro Institutionum calculi
integralis (publicado en 1768-1770).
El método de Euler es un método de primer orden, lo que significa que el error local
es proporcional al cuadrado del tamaño del paso, y el error global es proporcional al
tamaño del paso. El método de Euler regularmente sirve como base para construir
métodos más complejos.

Ilustración del método de Euler. La curva desconocida está en azul, y su

aproximación poligonal en rojo.

euler_20220425a22.m

Algoritmo:

ff=@(t,y) t-y;

y0=0; t=a; a=0; b=2; n=4;

h=(b-a)/n;

for i=1:n
y1=y0+h*ff(t,y0);
printf('%d t=%f y=%f\n',i,t+h,y1);
y0=y1;
t=t+h;

end;

%Aumentando Datos 1
 printf('-------------')

ff=@(t,y) t-y;

y0=0; t=a; a=0; b=2; n=4;

h=(b-a)/n; yt(1)=y0

for i=1:n
y1=y0+h*ff(t,y0); yt(i+1)=y1;
printf('%d t=%f y=%f\n',i,t+h,y1);
y0=y1;
t=t+h;
end;
tx=[0 1/2 1 3/2 2];
plot(tx,yt)

%Aumentando Datos 2
printf('-------------')

ff=@(t,y) t-y;

y0=0; t=a; a=0; b=2; n=4;

h=(b-a)/n; yt(1)=y0

for i=1:n
y1=y0+h*ff(t,y0); yt(i+1)=y1;
printf('%d t=%f y=%f\n',i,t+h,y1);
y0=y1;
t=t+h;
end;
tx=[0 1/2 1 3/2 2]; t=0:.01:2;y=t-1+exp(-t);
plot(tx,yt,t,y)

euler_20220425a55.m
Algoritmo:
ff=@(t,y) sqrt(t.^2+y.^(3/2));

y0=0; t=a; a=0; b=2; n=12;

h=(b-a)/n; yt(1)=y0

for i=1:n
y1=y0+h*ff(t,y0); yt(i+1)=y1;
printf('%d t=%f y=%f\n',i,t+h,y1);
y0=y1;
t=t+h;

3
 end;
tx=0:h:2;
plot(tx,yt)

euler_20220425compuesto11.m

Algoritmo:

ff=@(t,y) t-y;

y0=0; t=a; a=0; b=2; n=4;

h=(b-a)/n; yk=y0;
for i=1:n
for j=1:10
y1=y0+h/2*(ff(t,y0)+ff(t+h,yk));
if abs(y1-yk)<0.001
break;
end;
yk=y1;
end,
y0=y1; t=t+h;
printf('%d t=%f y=%f\n',i,t,y1);
end;

euler_20220425compuesto22.m

Algoritmo:

%arregos de compuesto11

ff=@(t,y) t-y;

y0=0; t=a; a=0; b=2; n=4;

h=(b-a)/n;
yt=[]; %ojo
yk=y0; yt(1)=y0;
for i=1:n
for j=1:10
y1=y0+h/2*(ff(t,y0)+ff(t+h,yk));
if abs(y1-yk)<0.001
break;
end;
yk=y1;
end,
y0=y1; t=t+h; yt(i+1)=y1;
printf('%d t=%f y=%f\n',i,t,y1);

4
end;
tx=0:h:2;
%para comparar
t=0:0.01:2; y=t-1+exp(-t); plot(tx,yt,t,y);

printf('------------')

ff=@(t,y) t-y;

y0=0; t=a; a=0; b=2; n=12;

h=(b-a)/n;
yt=[]; %ojo
yk=y0; yt(1)=y0;
for i=1:n
for j=1:10
y1=y0+h/2*(ff(t,y0)+ff(t+h,yk));
if abs(y1-yk)<0.001
break;
end;
yk=y1;
end,
y0=y1; t=t+h; yt(i+1)=y1;
printf('%d t=%f y=%f\n',i,t,y1);
end;
tx=0:h:2;
%para comparar
t=0:0.01:2; y=t-1+exp(-t); plot(tx,yt,t,y);

5
 MÉTODO DE ADAMS MOULTON

¿Qué es el Método de Método de Adams Moulton?

En este trabajo de investigación se construyó un algoritmo para solucionar un

problema de valor inicial con derivada fraccionária tipo Caputo. Para tal
construcción, se tomó como base el método de Adams-Bashforth-Moulton para
ecuaciones diferenciales ordinaria de orden entero. En el desarrollo se presenta tres
ejemplos numéricos, en los cuales, fue estimado numéricamente el orden de
convergencia, separando en dos grupos, el primero para el orden de la ecuación
diferencial fraccionária que varia entre 0 y 1 y el segundo que varia entre 1 y 2. En
cada grupo, se presenta una tabla con el cálculo explícito del orden de convergencia
para los diferentes valores del orden de la ecuación diferencial fraccionária. La
presente investigación es de tipo básica, y el método utilizado es el inductivo-
deductivo pues se analiza la teoría del método Adams-Bashforth-Moulton para una
ecuación diferencial ordinaria de primer orden y luego se hace una generalización de
dicho método para una ecuación diferencial fraccionária de tipo Caputo. En los
experimentos computacionales se consideran ciertas ecuaciones diferenciales de
orden fraccionário que se conocen explícitamente su solución y se comparan con la
solución obtenida con el método generalizado presentado. Al respecto se halla
intuitivamente que el orden de convergencia en ciertos casos particulares es
aproximadamente el minimo entre 2 y el orden de la ecuación diferencial más 1; y en
otro caso el orden de convergencia disminuye conforme aumenta el valor del orden
de la ecuación diferencial.

ecuacionesss.m

Algoritmo:

%20220516

f1=@(t,x,y) 2*x-3*y;
f2=@(t,x,y) x+y;
a=0; b=2; y0=1; x0=0; n=20; h=(b-a)/n;
t=a;
yp1(1)=x0;
%yp1=[];
yp2(1)=y0;
%yp2=[];
for i=1:n
k11=h*f1(t,x0,y0); k12=h*f2(t,x0,y0);
k21=h*f1(t+h/2, x0+k11/2, y0+k12/2); k22=h*f2(t+h/2, x0+k11/2, y0+k12/2);

6
 k31=h*f1(t+h/2, x0+k21/2, y0+k22/2); k32=h*f2(t+h/2, x0+k21/2, y0+k22/2);
k41=h*f1(t+h , x0+k31 , y0+k32) ; k42=h*f2(t+h , x0+k31 , y0+k32) ;
x1 = x0 + (k11 + 2*(k21+k31)+k41)/6;
y1 = y0 + (k12 + 2*(k22+k32)+k42)/6;

printf('%d t=%f x0=%f y0=%f\n',i , t+h, x1, y1);

yp1(i+1)=x1;
x0=x1;
yp2(i+1)=y1;
y0=y1;
t=t+h;
end;
xt=0:h:2;
yt=0:h:2;
%xt=[];
%yt=[];
plot(xt,yp1, 'r', 'linewidth', 1);
grid on
hold on
plot(yt,yp2, 'b', 'linewidth', 1);
grid on
hold on
plot(y=0, 's', 'linewidth', 1)
hold off

adamsa_Moulton_11.m

Algoritmo:

%20220516

% Adams- moulton

% n es el orden 2 o numero de pasos 2

% formula
% yn+1 = yn + (h/12)*(5*f(tn, yn+1)8*f(tn, yn)-f(tn-1, yn-1))
%
function adamsa_Moulton_11
clear all;
ff=@(t,y) t-y;
a=0; b=2; y0=0; n=8; h=(b-a)/n;
y1=runge_kutta_a1(a,1*h, y0, 1)
tt=2*h; yp=y1; y2=y1; yn1=y1(2);y1(3)=yn1;
for i=1:n-1
y1(3)=yn1;
for j=1:10 % circular 2 pasos
yn1=(5*ff(tt, y1(3))+8*ff(tt-h, y1(2))-ff(tt-2*h, y1(1)));
yn1=(yn1)*h/12+y1(2);
yr=abs(y1(3)-yn1);
y1(3)=yn1;
7
 if yr<.0001 break; end;
end;
%y1(3)=y2(4);
y1(2)=yn1;r22=y1(2)
y1(1)=y2(2); tt=tt+h;
y2=y1;yp(i+2)=yn1;
end;
x=0:h:2; t=0:.01:2; yw=exp(-t)+t-1;
plot(x,yp, t, yw);
grid;
end;%function moulton;
%%%%% funcion Runge Kutta

%runge_kutta_a
function yw=ff(t,y) yw=t-y; end;

function yt1=runge_kutta_a1(a,b, y0, n)

%ff=@(t,y) t-y;
%a=0; b=2; y0=0; n=4;
h=(b-a)/n; t=a; k=0; y(1)=y0;

for i=1:n
k1=h*ff(t,y0);
k2=h*ff(t+h/2, y0+k1/2);
k3=h*ff(t+h/2, y0+k2/2);
k4=h*ff(t+h, y0+k3);
y1=y0+(k1+2*(k2+k3)+k4)/6;
k=k+1; y(k+1)=y1; y0=y1;
t=t+h;
end;
yt1=y;
endfunction;

adamsa_Moulton_22.m

Algoritmo:

%MULTON DE 3 PASOS
%20220516
function adamsa_Moulton_22
clear all;
ff=@(t,y) t-y;
a=0; b=2; y0=0; n=4; h=(b-a)/n;
y1=runge_kutta_a1(a,2*h, y0, 2)
tt=3*h;
yp=y1;
y2=y1;
yn1=y1(3);
y1(4)=yn1;
for i=2:n-1
y1(4)=yn1;
for j=1:10 % circular 3 pasos
8
 yn1=(9*ff(tt, y1(4))+19*ff(tt-h, y1(3))-5*ff(tt-2*h, y1(2))+ff(tt-3*h,
y1(1)));
yn1=(yn1)*h/24+y1(3);
yr=abs(y1(4)-yn1);
y1(4)=yn1;
if yr<.0001 break; end;
end;
%y1(3)=y2(4);
y1(3)=yn1;
r22=y1(3)
y1(2)=y2(3);
y1(1)=y2(2);
tt=tt+h;
y2=y1;
yp(i+2)=yn1;
end;
x=0:h:2; t=0:.01:2; yw=exp(-t)+t-1;
plot(x,yp, t, yw);
%xlabel('x');
%ylabel('y');
%title('Gráfica Adams - bashfort de dos pasos');
grid;
end;%function moulton;
%%%%% funcion Runge Kutta

%runge_kutta_a
function yw=ff(t,y) yw=t-y; end;

function yt1=runge_kutta_a1(a,b, y0, n)

%ff=@(t,y) t-y;
%a=0; b=2; y0=0; n=4;
h=(b-a)/n; t=a; k=0; y(1)=y0;

for i=1:n
k1=h*ff(t,y0);
k2=h*ff(t+h/2, y0+k1/2);
k3=h*ff(t+h/2, y0+k2/2);
k4=h*ff(t+h, y0+k3);
y1=y0+(k1+2*(k2+k3)+k4)/6;
k=k+1; y(k+1)=y1; y0=y1;
t=t+h;
end;
yt1=y;
endfunction;

ecuaciones.m

Algoritmo:
%clc
f1=@(t,x,y) 2*x-3*y;
f2=@(t,x,y) x+y;
a=0; b=2; y0=1; x0=0; n=20; h=(b-a)/n;

9
t=a;
yp1(1)=x0;
yp2(1)=y0;
for i=1:n
k11=h*f1(t,x0,y0); k12=h*f2(t,x0,y0);
k21=h*f1(t+h/2, x0+k11/2, y0+k12/2); k22=h*f2(t+h/2, x0+k11/2, y0+k12/2);
k31=h*f1(t+h/2, x0+k21/2, y0+k22/2); k32=h*f2(t+h/2, x0+k21/2, y0+k22/2);
k41=h*f1(t+h , x0+k31 , y0+k32) ; k42=h*f2(t+h , x0+k31 , y0+k32) ;
x1 = x0 + (k11 + 2*(k21+k31)+k41)/6;
y1 = y0 + (k12 + 2*(k22+k32)+k42)/6;

printf('%2d t=%f x0=%f y0=%f\n',i , t+h, x1, y1);

yp1(i+1)=x1;
x0=x1;
yp2(i+1)=y1;
y0=y1;
t=t+h;
end;
tx=0:h:2;
ty=0:h:2;
%ty=[];
%tx=[];
plot(tx,yp1)
hold on
plot(ty,yp2)
hold off

adamsa_Moulton_tres.m

Algoritmo:

% Adams- moulton

% n es el orden 3 o numero de pasos 3

% formula
% yn+1 = yn + (h/24)*(9*f(tn, yn+1)+19*f(tn, yn)
% -5f(tn-1, yn-1) + f(tn-2, yn-2) )
%
function adamsa_Moulton_tres
clear all;
ff=@(t,y) t-y;
a=0; b=2; y0=0; n=4; h=(b-a)/n;
y1=runge_kutta_a1(a,2*h, y0, 2)
tt=3*h; yp=y1; y2=y1; yn1=y1(3);y1(4)=yn1;
for i=2:n-1
y1(4)=yn1;
for j=1:10 % circular 2 pasos
yn1=(9*ff(tt, y1(4))+19*ff(tt-h, y1(3))-5*ff(tt-2*h, y1(2)));

10
 yn1=(yn1+ff(tt-3*h, y1(1)))*h/24+y1(3);
yr=abs(y1(4)-yn1);
y1(4)=yn1;
if yr<.0001 break; end;
end;
%y1(3)=y2(4);
y1(3)=yn1;r22=y1(3)
y1(2)=y2(3); y1(1)=y2(2); tt=tt+h;
y2=y1;yp(i+2)=yn1;
end;
x=0:h:2; t=0:.01:2; yw=exp(-t)+t-1;
plot(x,yp, t, yw);
grid;
end;%function moulton;
%%%%% funcion Runge Kutta

%runge_kutta_a
function yw=ff(t,y) yw=t-y; end;

function yt1=runge_kutta_a1(a,b, y0, n)

h=(b-a)/n; t=a; k=0; y(1)=y0;
for i=1:n
k1=h*ff(t,y0);
k2=h*ff(t+h/2, y0+k1/2);
k3=h*ff(t+h/2, y0+k2/2);
k4=h*ff(t+h, y0+k3);
y1=y0+(k1+2*(k2+k3)+k4)/6;
k=k+1; y(k+1)=y1; y0=y1;
t=t+h;
end;
yt1=y;
endfunction;

11
12
 MÉTODO DE RUNGE KUTTA

Uno de los métodos más utilizados para resolver numéricamente problemas de

ecuaciones diferenciales ordinarias con condiciones iniciales es el método de Runge-
Kutta de cuarto orden, el cual proporciona un pequeño margen de error con respecto
a la solución real del problema y es fácilmente programable en un software para
realizar las iteraciones necesarias. El método de Runge-Kutta se utiliza para resolver
ecuaciones diferenciales de la forma:

Y es sumamente útil para casos en los que la solución no puede hallarse por los
métodos convencionales (como separación de variables). Hay variaciones en el
método de Runge-Kutta de cuarto orden pero el más utilizado es el método en el cual
se elige un tamaño de paso h y un número máximo de iteraciones n tal que

y se realiza la iteración Para i = 0,…, n-1. La solución se da a lo largo del intervalo (to,
to+ hn).

Runge_kuttaxx.m

Algoritmo:

13
function ww1=Runge_kuttaxx(a,b,ya, yp,n)
% Runge_kuttayy(a,b,ya_front, yp_penditne,ya_deriv, yp_deriv, n_pasos)
f1=@(t,x,y) 2*x-3*y; f2=@(t,x,y) x+y;
a=0; b=2; n=20; h=(b-a)/n; ya=0; yp=1;t=a;
h=(b-a)/n;
t=a;ww1=[];
ww1(1,1)=ya; ww1(2,1)=yp;
for ii=0:n-1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k11=h*f1(t,ya,yp);
k21=h*f2(t,ya,yp);
k12=h*f1(t+h/2,ya+1/2*k11,yp+1/2*k21);
k22=h*f2(t+h/2,ya+1/2*k11,yp+1/2*k21);
k13=h*f1(t+h/2,ya+1/2*k12,yp+1/2*k22);
k23=h*f2(t+h/2,ya+1/2*k12,yp+1/2*k22);
k14=h*f1(t+h,ya+k13,yp+k23);
k24=h*f2(t+h,ya+k13,yp+k23);
ya1=ya+(k11+2*k12+2*k13+k14)/6;
y11=yp+(k21+2*k22+2*k23+k24)/6;
ww1(1,ii+2)=ya1;ww1(2,ii+2)=y11;
t=t+h;ya=ya1; yp=y11;
end;
[ww1']
t=0:h:2;
plot(t, ww1(1,:), t, ww1(2,:));

runge_kutta_a.m

Algoritmo:

%runge_kutta_a

function yt1=runge_kutta_a(a,b, y0, n)

ff=@(t,y) t-y;
%a=0; b=2; y0=0; n=4;
h=(b-a)/n; t=a; k=0; y(1)=y0;

for i=1:n
k1=h*ff(t,y0);
k2=h*ff(t+h/2, y0+k1/2);
k3=h*ff(t+h/2, y0+k2/2);
k4=h*ff(t+h, y0+k3);
y1=y0+(k1+2*(k2+k3)+k4)/6;
k=k+1; y(k+1)=y1; y0=y1;
t=t+h;
end;
yt1=y;

14
15
 MÉTODO DE DIFERENCIAS FINITAS

¿Qué es el Método de Diferencias Finitas?

Una diferencia finita es una expresión matemática de la forma f(x + b) − f(x +a). Si
una diferencia finita se divide por b − a se obtiene una expresión similar
al cociente diferencial, que difiere en que se emplean cantidades finitas en lugar
de infinitesimales. La aproximación de las derivadas por diferencias finitas
desempeña un papel central en los métodos de diferencias finitas del análisis
numérico para la resolución de ecuaciones diferenciales.

16
 diferencias_finitas_11_20220602.m

Algoritmo:

A=[-(2+(10/4)^2 * 1/2) 1 0; 1 -(2+(10/4)^2 * 1/2) 1; 0 1 -

(2+(10/4)^2 * 1/2)];

B=[-(10/4)^2*exp(-2.5/5)/2-1; -(10/4)^2*exp(-5/5)/2; -
(10/4)^2*exp(-7.5/5)/2-4];

c=A\B

Difere_finita1.m

Algoritmo:

% EL VALOR DE k SI ES CERO LE PONE LOS PASOS INTEMEDIOS SI ES 1 NO

function d=Difere_finita1(a,b,ya,yb,n,k)
d=[];
AA=[];
a=1; b=3; ya=17; yb=43/3;n=20; k=2;
h=(b-a)/n;
AA(1,1)=2+h*h*qq(a+h);
AA(1,2)=-1+h*pp(a+h)/2;
RR(1)=-h*h*rr(a+h)+(1+h*pp(a+h)/2)*ya;
for i=2:n-2
AA(i,i-1)=-1-h*pp(a+h*i)/2;
AA(i,i)=2+h*h*qq(a+h*i);
AA(i,i+1)=-1+h*pp(a+h*i)/2;
RR(i)=-h*h*rr(a+i*h);
end
AA(n-1,n-2)=-1-h*pp(a+h*(n-1))/2;
AA(n-1,n-1)=2+h*h*qq(a+h*(n-1));
17
RR(n-1)=-h*h*rr(a+h*(n-1))+(1-h*pp(a+h*(n-1))/2)*yb;
if k<1
AA
RR
end
format long
s=AA\RR'
format;
plot(s); grid;
endfunction;

function t1=pp(x)
%t1=0;
t1=-2/x;
end;
function t2=qq(x)
%t2=1/2;
t2=2/x^2;
end;
function t3=rr(x)
%t3=-exp(-0.2*x)/2;
t3=sin( log(x) ) / x^2;
end;

Difere_finita33.m

Algoritmo:

% EL VALOR DE k SI ES CERO LE PONE LOS PASOS INTEMEDIOS SI ES 1 NO

% y" = p(x)y' + q(x)y + r(x)
% y(a)=ya; y(b) = yb; a<= x <= b;
function d=Difere_finita33(a,b,ya,yb,n,k)
d=[];
AA=[];
%valores a b y ya y yb mas n=pasos
a=0; b=2; n=10; ya=0; yb=1; k=.2;
h=(b-a)/n;
AA(1,1)=2+h*h*qq(a+h);
AA(1,2)=-1+h*pp(a+h)/2;
RR(1)=-h*h*rr(a+h)+(1+h*pp(a+h)/2)*ya;
for i=2:n-2
AA(i,i-1)=-1-h*pp(a+h*i)/2;
AA(i,i)=2+h*h*qq(a+h*i);
AA(i,i+1)=-1+h*pp(a+h*i)/2;
RR(i)=-h*h*rr(a+i*h);
end
AA(n-1,n-2)=-1-h*pp(a+h*(n-1))/2;
AA(n-1,n-1)=2+h*h*qq(a+h*(n-1));
RR(n-1)=-h*h*rr(a+h*(n-1))+(1-h*pp(a+h*(n-1))/2)*yb;
if k<1
AA
RR
18
end
format long
s=AA\RR'; s';
format;
x=a:h:b; s1=[ya s' yb];
plot(x, s1, 'linewidth',3); grid;

end; % fin Difere_finita22

% aqui funciones p(x), q(x), r(x)
function t1=pp(x)
t1=-y/2
%t1=0;
%t1=-2/x;
end;
function t2=qq(x)
t2=0
%t2=1/10;
% t2=2/x^2;
end;
function t3=rr(x)
t3=x/2
%t3=x
%t3=-exp(-0.2*x)/2;
%t3=sin( log(x) ) / x^2;
end;

Difere_finita22.m

Algoritmo:

% EL VALOR DE k SI ES CERO LE PONE LOS PASOS INTEMEDIOS SI ES 1 NO

% y" = p(x)y' + q(x)y + r(x)
% y(a)=ya; y(b) = yb; a<= x <= b;
function d=Difere_finita22(a,b,ya,yb,n,k)
d=[];
AA=[];
%valores a b y ya y yb mas n=pasos
a=1; b=2; n=3; ya=1.5; yb=3; k=.5;
h=(b-a)/n;
AA(1,1)=2+h*h*qq(a+h);
AA(1,2)=-1+h*pp(a+h)/2;
RR(1)=-h*h*rr(a+h)+(1+h*pp(a+h)/2)*ya;
for i=2:n-2
AA(i,i-1)=-1-h*pp(a+h*i)/2;
AA(i,i)=2+h*h*qq(a+h*i);
AA(i,i+1)=-1+h*pp(a+h*i)/2;
RR(i)=-h*h*rr(a+i*h);
end
AA(n-1,n-2)=-1-h*pp(a+h*(n-1))/2;
AA(n-1,n-1)=2+h*h*qq(a+h*(n-1));
RR(n-1)=-h*h*rr(a+h*(n-1))+(1-h*pp(a+h*(n-1))/2)*yb;
if k<1
19
 AA
RR
end
format long
s=AA\RR'; s';
format;
x=a:h:b; s1=[ya s' yb];
plot(x, s1, 'linewidth',2); grid;

end; % fin Difere_finita22

% aqui funciones p(x), q(x), r(x)
function t1=pp(x)
t1=-1/x;
%t1=0;
%t1=-2/x;
end;
function t2=qq(x)
t2=0;
%t2=1/2;
% t2=2/x^2;

end;
function t3=rr(x)
t3=0;
%t3=-exp(-0.2*x)/2;
%t3=sin( log(x) ) / x^2;
end;

20
21
 MÉTODO DEL DISPARO LINEAL

¿Qué es el Método del Disparo Lineal?

Una EDO lineal con condiciones de borde es de la forma:

El método del disparo para EDO lineales, consiste en sustuir el problema lineal con
valor de frontera por dos problemas de valor inicial. Por lo tanto se resolverán los
siguientes problemas de valor inicial:

La solución puede construirse entonces usando las

soluciones e como:

disparo_lineal_11.m
Algoritmo:

function yy1=disparo_lineal_11
a=0; b=1; n=2; h=(b-a)/n;
yA=0;yB=1;

YY1=Runge_kuttaxx(a, b, yA,0, n);

YY2=Runge_kuttaxx(a, b, 0, n);

YY1(1, n+1), YY2(1, n+1)

YY=YY1(1,:) + (yB-YY1(1, n+1))/YY2(1, n+1).*YY2(1,:)
end;

function y=f1(t y0, w);

y=w:
end;

22
 MÉTODO DEL DISPARO NO LINEAL

¿Qué es el Método del Disparo no Lineal?

Para las EDOs no lineales, no es posible expresar la solución como una

combinación lineal de las soluciones de dos problemas de valor inicial, sino que se
necesita un conjunto de soluciones que permitan ir convergiendo a la solución.
El método del disparo permite resolver EDOs no lineales de segundo orden de la
forma:

que consiste en resolver el problema de valor inicial:

(1)
para varios valores del parámetros hasta obtener una solución que
satisfaga .
Primero definamos la función:

23
 notar que esta función es positiva cuando se está sobreestimando el borde y es
negativa cuando el valor de es menor que .
Cuando encontramos dos valores para digamos y que cumplen lo
siguiente:

entonces podemos decir que podemos encontrar un valor tal que que
está entre e .

La manera de estimar a partir de y es usando el método de regula falsi:

disparo_no_lineal_secante11_20220530.m

Algoritmo:

% DISPARO NO LINEAL

function y1 = disparo_no_lineal_secante11_20220530
a=0; b=2; ya=0; wa=1; n=10; p=1; p1=0; yb=1; % varia

ys = rungekuttageneral(a, b, ya, p, n); %primera llamada

ys1 = ys(length(ys));
for i=1:80
ys = rungekuttageneral(a, b, ya, p1, n);
ys2 = ys(length(ys));
% secante
p2=p-(ys1-1)*(p-p1)/(ys1-ys2);
disp([p2 ys2]);
p4=abs(p1-p);
p1=p; p=p2; ys1=ys2;
24
 if p4<0.00001 break; end;
end;
ys
plot(ys);
grid;
end; %fin de disparo ...

function y2=f1(t, y, w)
y2=w;
end;
function y2=f2(t, y, w)
y2=t+y*y/10;
end;
%runge
function yt1 = rungekuttageneral(a, b, ya, wa, n)
h = (b-a)/n;
t=a; k=0; %y(1)=y0;

for i=1:n
ky1 = h*f1(t,ya,wa);
kw1 = h*f2(t,ya,wa);
ky2 = h*f1(t+h/2,ya+ky1/2, wa+kw1/2);
kw2 = h*f2(t+h/2,ya+ky1/2, wa+kw1/2);
ky3 = h*f1(t+h/2,ya+ky2/2, wa+kw2/2);
kw3 = h*f2(t+h/2,ya+ky2/2, wa+kw2/2);
ky4 = h*f1(t+h,ya + ky3, wa+kw3);
kw4 = h*f2(t+h,ya + ky3, wa+kw3);

y1=ya+(ky1 + 2*ky2 + 2*ky3 + ky4)/6;

y2=wa+(kw1 + 2*kw2 + 2*kw3 + kw4)/6;

k=k+1; y(k+1)=y1; ya=y1; wa=y2;

t = t+h;
end;
yt1=y;
end;

25
 MÉTODO DE RAYLEIGH-RITZ

¿Qué es el Método de Rayleigh-Ritz?

El método de Rayleigh-Ritz es un método numérico directo de aproximación

de valores propios , originado en el contexto de la resolución de problemas físicos
de valores límite y llamado así por Lord Rayleigh y Walther Ritz .
Según, citando a Richard Courant , tantoLord Rayleigh y Walther Ritz concibieron de
forma independiente la idea de utilizar la equivalencia entre los
problemas de valores en la frontera de las ecuaciones diferenciales parciales, por un
lado, y los problemas de cálculo de variaciones , por otro lado, para el cálculo
numérico de las soluciones, sustituyendo los problemas variacionales más simples.
aproximar problemas extremos en los que es necesario determinar un número finito
de parámetros; consulte el artículo Método Ritz para obtener más
información. Irónicamente para el debate, la justificación moderna del algoritmo
abandona el cálculo de variaciones a favor del enfoque más simple y general
de proyección ortogonal como en el método Galerkin lleva el nombre de Boris
Galerkin , lo que lleva también a la nomenclatura del método Ritz-Galerkin .

26
 Se utiliza en todas las aplicaciones que implican la aproximación de valores propios y
vectores propios , a menudo con nombres diferentes. En la mecánica cuántica , donde
se describe un sistema de partículas mediante un hamiltoniano , el método de
Ritz utiliza funciones de onda de prueba para aproximar la función propia del estado
fundamental con la energía más baja. En el contexto del método de elementos finitos ,
matemáticamente el mismo algoritmo se denomina comúnmente método Ritz-
Galerkin . El método Rayleigh-Ritz o la terminología del método Ritz es típica en
ingeniería mecánica y estructural para aproximar los modos propios y las frecuencias
resonantes.de una estructura.

Rayleigh_20220630_11_frontera__.m

Algoritmo:

% programa hecho por "andres" y corregido

% -(p(x)*y')'+q(x)*y=f(x)
% y''-4y=x ; -(1*y')'+4y=-x;

% TRABAJA PARA [-2, 3] Y LO GRAFICA

% VEA LOS PORMENORES EN EL VIDEO
% PARA QUE ENTIENDA
# OK
clear
p=@(x)1; q=@(x)4*25; f=@(x)-25*(5*x-2);
x=[0 .05 .1 .2 1/3 .4 .5 .6 2/3 .7 .8 .9 .95 1]; n=length(x)- 2;

h=x(2)-x(1);
Ip1=@(x)p(x); Ip1=quad(Ip1,x(1),x(2))/h^2;

Iq1=@(z)(z-x(1)).^2*q(z); Iq1=quad(Iq1,x(1),x(2))/h^2;

for i=1:n-1

h=x(i+2)-x(i+1);
IpN=@(z)p(z); IpN=quad(IpN,x(i+1),x(i+2))/h^2;
IqN=@(z)(z-x(i+2)).^2*q(z); IqN=quad(IqN,x(i+1),x(i+2))/h^2;
Ia=@(z)(x(i+2)-z).^2*q(z); Ia= quad(Ia,x(i+1),x(i+2))/h^2;
Ib=@(z)(x(i+2)-z).*(z-x(i+1)).*q(z); Ib= quad(Ib,x(i+1),x(i+2))/h^2;

h1=h;

h=x(i+1)-x(i);
Ic=@(z)(z-x(i)).*f(z); Ic= quad(Ic,x(i),x(i+1))/h;
Id=@(z)(x(i+2)-z).*f(z); Id= quad(Id,x(i+1),x(i+2))/h1;

A(i,i)=Ip1+IpN+Iq1+Ia;
A(i,i+1)=-IpN+Ib;
A(i+1,i)=A(i,i+1);
27
b(i)=Ic+Id;

Ip1=IpN; Iq1=IqN;
endfor
i=n;
h=x(i+2)-x(i+1);
IpN=@(z)p(z); IpN=quad(IpN,x(i+1),x(i+2))/h^2;
%IqN=@(z)(z-x(i+2)).^2*q(z); IqN=quad(IqN,x(i+1),x(i+2))/h^2;
Ia=@(z)(x(i+2)-z).^2*q(z); Ia= quad(Ia,x(i+1),x(i+2))/h^2;
Ib=@(z)(x(i+2)-z).*(z-x(i+1)).*q(z); Ib= quad(Ib,x(i+1),x(i+2))/h^2;

h1=h;

h=x(i+1)-x(i);
Ic=@(z)(z-x(i)).*f(z); Ic= quad(Ic,x(i),x(i+1))/h;
Id=@(z)(x(i+2)-z).*f(z); Id= quad(Id,x(i+1),x(i+2))/h1;
A(i,i)=Ip1+IpN+Iq1+Ia;
b(i)=Ic+Id;
A
b

c=A\b'
c=[0 c' 0];
% prporcion de la x =5*x-2;
x =5*x-2;
t=-2:.01:3;
%y=sinh(2*t)/(4*sinh(2))-t/4;
%c2=(2-exp(2))/(4*(exp(-4)-1)); c1=1/4*e^-2-c2*e^-4;
c1=(2*e^-10+3)/(4*(e^6-e^-14));
c2=-2/4*e^-4 -c1*e^-8;
y=c1*e.^(2*t)+c2*e.^(-2*t)-t/4;
plot(x, c, t, y); grid;

MetodoRayleighRitz23062022.m

Algoritmo:

%y''-4y=x; -(1y')'+4y=-x;
p=@(x) 1; q=@(x) 4; r=@(x)-x;
x=[0, 1/3, 2/3, 1]; n=2;
Ip1=@(z) p(z); h=x(2)- x(1); Ip1=quad(Ip1,x(1),x(2))/h^2
Iq1=@(z) (z-x(1)).^2.*q(z); Iq1=quad(Iq1,x(1),x(2))/h^2
for i=1:n-1
IpN=@(z) p(z);h=x(i+2)-x(i+1);IpN=quad(IpN,x(i+1),x(i+2))/h^2;
IqN=@(z) (z-x(i+2)).^2.*q(z);IqN=quad(IqN,x(i+1),x(i+2))/h^2;
Ia=@(z) (x(i+2)-z).^2.*q(z);Ia=quad(Ia,x(i+1),x(i+2))/h^2;
Ib=@(z) (x(i+2)-z).*(z-x(i+1));Ib=quad(Ib,x(i+1),x(i+2))/h^2;h1=h;
Ic=@(z) (z-x(i)).*r(z);h=x(i+1)-x(i);Ic=quad(Ic,x(i),x(i+1))/h;
Id=@(z) (x(i+2)-z).*r(z);Id=quad(Id,x(i+1),x(i+2))/h1;
A(i,i)=Ip1+IpN+Iq1+Ia; A(i,i+1)=-IpN+Ib; A(i+1,i)=A(i,i+1);
28
 b(i)=Ic+Id; Ip1=IpN;Iq1=IqN;
end;
i=n+1;
IpN=@(z) p(z);h=x(i)-x(i+1);IpN=quad(IpN,x(i-1),x(i))/h^2;
Ia=@(z) (x(i)-z).^2.*q(z);Ia=quad(Ia,x(i-1),x(i))/h^2;
Ib=@(z) (x(i)-z).*(z-x(i-1)).*q(z);Ib=quad(Ib,x(i-1),x(i))/h^2;
Ic=@(z) (z-x(i-1)).*r(z);h1=x(i-1)-x(i-2);Ic=quad(Ic,x(i-2),x(i-1))/h;
Id=@(z) (x(i)-z).*r(z);Id=quad(Id,x(i-1),x(i))/h;
A(n,n)=Ip1+IpN+Iq1+Ia;
b(n)=Ic+Id
A
b
c=A\b'

raleygth22_1.m

Algoritmo:

% y" - 2y = -x; y(0)=0; y(1)=0;

% y=-sinh(2^.5*x)/sinh(2^.5)/2 +x/2;
function ytw=raleygth22_1

x=[0 1/3 1/2 .6 2/3 .82 1];

x=0:1/7:1;x=[0 1/3 2/3 1];
%x=[0 .15 .3 .4 .45 .5 .55 .6 .65 .7 .8 .85 .9 1];
%x=[0 13/100 1/4 2/5 1/1.9 2/3 .76 .93 1];
pp='1'; qq='2'; ff='x';
n=length(x)-1;

Ip1s=[pp];Ip1f=inline(Ip1s);
Ip1=quad(Ip1f, x(1),x(2)).*(b1=1/(x(2)-x(1))^2);
Iq1s=['(' qq ').*(x-' mat2str(x(1)) ').^2'];Iq1f=inline(Iq1s);
Iq1=quad(Iq1f, x(1),x(2)).*b1;

for j=1:n-2 i=j+1; j;

%IpNs=['1'];IpNf=inline(IpNs);
IpN=quad(Ip1f, x(i),x(i+1)).*(b1=1/(x(i+1)-x(i))^2);
IqNs=['(' qq ').*(x-' mat2str(x(i)) ').^2'];IqNf=inline(IqNs);
IqN=quad(IqNf, x(i),x(i+1)).*b1;
Ias=['(' qq ').*(' mat2str(x(i)) '-x).^2'];Iaf=inline(Ias);
Ia=quad(Iaf, x(i),x(i+1)).*b1;
Ibs=['(' qq ').*(' mat2str(x(i+1)) '-x).*(x-' mat2str(x(i))
')'];Ibf=inline(Ibs);
Ib=quad(Ibf, x(i),x(i+1)).*b1;
Ics=['(' ff ').*(x-' mat2str(x(i-1)) ')'];Icf=inline(Ics);
Ic=quad(Icf, x(i-1),x(i))/(x(i)-x(i-1));
Ids=['(' ff ').*(' mat2str(x(i+1)) '-x)'];Idf=inline(Ids);
Id=quad(Idf, x(i),x(i+1))/(x(i+1)-x(i));
A(j,j)=Ip1+IpN+Iq1+Ia;
A(j,j+1)=-IpN+Ib;
A(j+1,j)=A(j,j+1);
b(j)=Ic+Id;
29
 Ip1=IpN; Iq1=IqN;

end;
i=n;
IpN=quad(Ip1f, x(i),x(i+1)).*(b1=1/(x(i+1)-x(i))^2);
Ias=['(' qq ').*(' mat2str(x(i+1)) '-x).^2'];Iaf=inline(Ias);
Ia=quad(Iaf, x(i),x(i+1)).*b1;
%Ibs=['2*(' mat2str(x(i+1)) '-x).*(x-' mat2str(x(i)) ')'];Ibf=inline(Ibs);
%Ib=quad(Ibf, x(i),x(i+1)).*b1
Ics=['(' ff ').*(x-' mat2str(x(i-1)) ')'];Icf=inline(Ics);
Ic=quad(Icf, x(i-1),x(i))/(x(i)-x(i-1));
Ids=['(' ff ').*(' mat2str(x(i+1)) '-x)'];Idf=inline(Ids);
Id=quad(Idf, x(i),x(i+1))/(x(i+1)-x(i));
A(i-1,i-1)=Ip1+IpN+Iq1+Ia
b(i-1)=Ic+Id
c=A\b'; c'
y=[0 c' 0];
t=0:.01:1; yt=-sinh(2^.5*t)/sinh(2.^.5)/2 +t/2;
plot(t,yt,'linewidth',2, x, y,'linewidth',2);
grid;

rayleigh-ritz.m

Algoritmo:

% -(p(x)*y')'+q(x)*y=f(x)
% y''-4y=x ; -(1*y')'+4y=-x;
clear
p=@(x)1; q=@(x)4; f=@(x)(-x);
x=[0 1/3 2/3 1]; n=2;

h=x(2)-x(1);
Ip1=@(x)p(x); Ip1=quad(Ip1,x(1),x(2))/h^2;

Iq1=@(z)(z-x(1)).^2*q(z); Iq1=quad(Iq1,x(1),x(2))/h^2;

for i=1:n-1

h=x(i+2)-x(i+1);
IpN=@(z)p(z); IpN=quad(IpN,x(i+1),x(i+2))/h^2;
IqN=@(z)(z-x(i+2)).^2*q(z); IqN=quad(IqN,x(i+1),x(i+2))/h^2;
Ia=@(z)(x(i+2)-z).^2*q(z); Ia= quad(Ia,x(i+1),x(i+2))/h^2;
Ib=@(z)(x(i+2)-z).*(z-x(i+1)); Ib= quad(Ib,x(i+1),x(i+2))/h^2;

h1=h;

h=x(i+1)-x(i);
Ic=@(z)(z-x(i)).*f(z); Ic= quad(Ic,x(i),x(i+1))/h;
Id=@(z)(x(i+2)-z).*f(z); Id= quad(Id,x(i+1),x(i+2))/h1;

A(i,i)=Ip1+IpN+Iq1+Ia;
A(i,i+1)=-IpN+Ib;
30
A(i+1,i)=A(i,i+1);

b(i)=Ic+Id;

Ip1=IpN; Iq1=IqN;
endfor
i=n;

31
32
33
 MÉTODO DE FIBONACCI

 ¿Qué es el Método de Fibonacci?

Se trata de una secuencia infinita de números naturales; a partir del 0 y el 1, se van

sumando a pares, de manera que cada número es igual a la suma de sus dos
anteriores, de manera que:
0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55…

metodoFibonacci.m

Algoritmo:

%%%%%%%%%%%%%%%%%%%%%% método de fibonacci %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% para el problema multidimensional

function [alpha1]=metodoFibonacci(f,x0,Fgrad,bk,Max)

ak=0; % valor inicial

tol=0.03; % tolerancia
%% definición del usuario
bk=5; % límite superior el usuario debe definirlo

%% inicio del algoritmo

N=(bk-ak)/tol; % se define el número de iteraciones basadas
en la tolerancia
Fib=[1 1]; % datos de inicio de la secuencia
i=2; % i=2 para empezar a llenar el vector de
fibonacci
while N>=Fib(end) % cuando se halla superado el valor de N
Fib(i+1)=Fib(i)+Fib(i-1); % serie de fibonacci
i=i+1; % se incrementa el contador
end

lamda=ak+(Fib(end-2)/Fib(end))*(bk-ak); % lamda
miu=ak+(Fib(end-1)/Fib(end))*(bk-ak); % miu
fl=f(x0(1,1)+lamda*Fgrad(1,1),x0(2,1)+lamda*Fgrad(2,1));
% se evalua lamda en la función
fm=f(x0(1,1)+miu*Fgrad(1,1),x0(2,1)+miu*Fgrad(2,1)); %
se evalua miu en la función
iter=0;
while size(Fib,2)>2

iter=iter+1; % incrementar el contador de las iteraciones

34
 if Max==0
if fl>fm % si la función en lamda es mayor que la
función de miu
ak=lamda; % actualizar el límite ak con lamda
lamda=miu; % actualizar lamda con el valor antiguo de
miu
fl=fm; % actualizar la evaluación de fl con el
valor de fm
miu=ak+(Fib(end-1)/Fib(end))*(bk-ak); % definir el nuevo valor de miu
fm=f(x0(1,1)+miu*Fgrad(1,1),x0(2,1)+miu*Fgrad(2,1));
% se evalua miu en la función
elseif fl<fm % si la función en miu es mayor que la
función de lamda
bk=miu % actualizar el límite bk con miu
miu=lamda; % actualizar miu con el valor antiguo de lamda
fm=fl; % actualizar la evaluación de fm con el
valor de fl
lamda=ak+(Fib(end-2)/Fib(end))*(bk-ak); % definir el nuevo valor de lamda
fl=f(x0(1,1)+lamda*Fgrad(1,1),x0(2,1)+lamda*Fgrad(2,1));
% se evalua lamda en la función
end
end
if Max==1
if fl>fm % si la función en lamda es mayor que la
función de miu
bk=miu % actualizar el límite bk con miu
miu=lamda; % actualizar miu con el valor antiguo de lamda
fm=fl; % actualizar la evaluación de fm con el
valor de fl
lamda=ak+(Fib(end-2)/Fib(end))*(bk-ak); % definir el nuevo valor de lamda
fl=f(x0(1,1)+lamda*Fgrad(1,1),x0(2,1)+lamda*Fgrad(2,1));
% se evalua lamda en la función
elseif fl<fm % si la función en miu es mayor que la
función de lamda
ak=lamda; % actualizar el límite ak con lamda
lamda=miu; % actualizar lamda con el valor antiguo de
miu
fl=fm; % actualizar la evaluación de fl con el
valor de fm
miu=ak+(Fib(end-1)/Fib(end))*(bk-ak); % definir el nuevo valor de miu
fm=f(x0(1,1)+miu*Fgrad(1,1),x0(2,1)+miu*Fgrad(2,1));
% se evalua miu en la función
end
end
Fib(end)=[];

D(iter,:)=[ak lamda miu bk ]; % registro de los intervalos

plot(lamda,fl,'*b',miu,fm,'*r')
hold on

end

35
alpha1=(lamda+miu)/2;

MÉTODO DE NEWTON MULTIDIMENSIONAL

¿Qué es el Método de Newton Multidimensional?

El método consiste en derivar parcialmente la función con respecto a las variables

presentes y después igualar a la función original multiplicada por un menos, y
repetir el procedimiento con todas las funciones del sistema. Los cuales son los
nuevos valores iniciales para X y Y.

NewtonMultidimensional.m

Algoritmo:

%%%%%%%%%%%%%%%%%%% Newton multidimensional %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
clc
%% datos que deben ser agregados por el usuario

% f=@(x1,x2) 2*x1*x2+2*x2-x1^2-2*x2^2; % función multidimensional

% grad={@(x1,x2) 2*x2-2*x1, % gradiente de la función
% @(x1,x2) 2*x1+2-4*x2};
% Hess={@(x1,x2) -2,@(x1,x2) 2 % gradiente de la función
% @(x1,x2) 2,@(x1,x2) -4};

f=@(x1,x2) (x1-2)^4+(x1-2*x2)^2; % función multidimensional

grad={@(x1,x2) 4*(x1-2)^3+2*(x1-2*x2), % gradiente de la función

@(x1,x2) -4*x1+8*x2};

Hess={@(x1,x2) 12*(x1-2)^2+2, @(x1,x2) -4 % hessiana de la función

@(x1,x2) -4, @(x1,x2) 8};

x0=[0;3]; % punto inicial debe ser vector columna

Tol=0.09; % tolerancia del algoritmo
MaxIter=100; % número de iteraciones máximas
% bk=5; % límite superior para la búsqueda
unidimensional

%% Inicio del algoritmo

Fgrad=[grad{1}(x0(1,1),x0(2,1)) % se evalua el gradiente en el punto

inicial

36
 grad{2}(x0(1,1),x0(2,1))];

Fhess=[Hess{1,1}(x0(1,1),x0(2,1)) Hess{1,2}(x0(1,1),x0(2,1));
Hess{2,1}(x0(1,1),x0(2,1)) Hess{2,2}(x0(1,1),x0(2,1))];

iter=1;
D(iter,:)={iter,x0,Fgrad,Fhess}
while norm(Fgrad)>=Tol && iter<=MaxIter

x0=x0-inv(Fhess)*Fgrad; % se calcula el nuevo punto

Fgrad=[grad{1}(x0(1,1),x0(2,1)) % se evalua el gradiente en el

punto
grad{2}(x0(1,1),x0(2,1))];

Fhess=[Hess{1,1}(x0(1,1),x0(2,1)) Hess{1,2}(x0(1,1),x0(2,1)); % se evalua la

hessiana en el punto
Hess{2,1}(x0(1,1),x0(2,1)) Hess{2,2}(x0(1,1),x0(2,1))];

iter=iter+1;
D(iter,:)={iter,x0,Fgrad,Fhess}
end

x0=x0-inv(Fhess)*Fgrad; % se calcula el nuevo punto

Fgrad=[grad{1}(x0(1,1),x0(2,1)) % se evalua el gradiente en el
punto
grad{2}(x0(1,1),x0(2,1))];

Fhess=[Hess{1,1}(x0(1,1),x0(2,1)) Hess{1,2}(x0(1,1),x0(2,1)); % se evalua la

hessiana en el punto
Hess{2,1}(x0(1,1),x0(2,1)) Hess{2,2}(x0(1,1),x0(2,1))];
iter=iter+1;
D(iter,:)={iter,x0,Fgrad,Fhess}

% f=[1 -1]
% A=[-16 8]
% B=25;
% lb=[-2;-2]
% ub=[2;2]
% linprog(f,A,B,[],[],lb,ub)

base_segunco11.m

Algoritmo:

% programa para {1 x x^2 x^3 . . . }

%EDO y'' -4y = x; y(0)=0; x=1/3 2/3

function ya=base_segunco11
a=0; b=1; n=5; h=(b-a)/n;

37
x1=a:h:b
A=[]; b=[];
% construcion de la matriz A
for i=1:n+1
A(1,i)=x_x(x1(1),i-1);
A(n+1,i)=dx_x(x1(1),i-1);
b(i)=rr(x1(i));
end;
for i=2:n
for j=1:n+1
k = ddx_x( x1(i),j-1);

A(i,j)=k+0*dx_x(x1(i),j-1)-4*x_x(x1(i),j-1);
end;
end;
A
b
% cosntruimos b esta arriba
format rat;
c=A\b'
format;
x=0:.1:1; ys=[];

%ys=c(1)+c(2)*x+c(3)*x.^2+c(4)*x.^3;
for i=1:length(x)
k=0;
for j=1:length(c)
k=k+c(j)*x(i).^(j-1);
endfor;
ys(i)=k;

end;

yt=5*sinh(2*x)/8-x/4;

plot(x, ys, x, yt);

grid;

% funciones direxct, deriv, deriv_deriv

function v=x_x(t,m)
if m==0; v=1; else v=t.^m; end;

function v=dx_x(t,m)
if m<=0 v=0;
elseif m==1 v=1;
else v=m*t.^(m-1);

38
 endif
function v=ddx_x(t,m)
if m<=1 v=0;
elseif m==2 v=2;
else v=m*(m-1)*t.^(m-2);
end;
function v=rr(t)
v=t;

ECUACIÓN HIPERBÓLICA PARA UNA EDP

¿Qué es la Ecuación Hiperbólica para una EDP?

Una ecuación hiperbólica en derivadas parciales es una ecuación diferencial en

derivadas parciales de segundo orden del tipo:

en la cual la matriz:

cuyos coeficientes pueden ser constantes o funciones

continuas en las variables (x,y), tiene
un determinante negativo.
Un ejemplo de una ecuación diferencial en derivadas parciales hiperbólica es
la ecuación de ondas:

hiperbolica11.m

Algoritmo:

% Ecuacion HIPRBOLICA
%
%
% necesitamos crear el vector inicial w0={w(1,0) w(2,0) w(3,0) ... w(n-1,0) }
% tambien Ut(i,0) i=1 2 ... n-1
% w (i,0) =f(xi ), i=1 2 ... n-1
%
% w(i,1) = u(i,0) + kkg(xi), i=1 2 ... n-1
%
% para matriz A = [ λ*λ*w(i?ˆ’1,j) + 2(1?ˆ’λ*λ)w(i,j) + λ*λ*w(i+1,j)]
% λ = a*kk/hh ;

39
% w(i,j+1) = 2(1?ˆ’λ*λ)w(i,j) + λ*λ*( w(i+1,j)+w(i?ˆ’1,j))-w(i,j-1)
% i=1 2 ... n-1
%
clear all;
tt=1;
n=10; % numero de psrticiones para x
a=0; b=1;
hh=(b-a)/n; % el h para x
kk=.01; % Ojo es para llegar a t>0; t=1.0.
m=tt/kk;
if floor(m)<m m=floor(m+1); end;
kk=tt/m;
% necesitamos crear el vector inicial w0={w(1,0) w(2,0) w(3,0) ... w(n-1,0) }
% tambien Ut(i,0) i=1 2 ... n-1
% w (i,0) =f(xi ), i=1 2 ... n-1
%
% w(i,1) = u(i,0) + kkg(xi), i=1 2 ... n-1
%
% para matriz A = [ λ*λ*w(i?ˆ’1,j) + 2(1?ˆ’λ*λ)w(i,j) + λ*λ*w(i+1,j)]
% λ = a*kk/hh ;
% w(i,j+1) = 2(1?ˆ’λ*λ)w(i,j) + λ*λ*( w(i+1,j)+w(i?ˆ’1,j))-w(i,j-1)
% i=1 2 ... n-1
aa=4;
x=a+hh:hh:b-hh;
ff=inline('sin(pi*x)'); % condicion inicial
w00=ff(x);
gg=inline('x*0'); % codicion de primera derivada de U
w11=w00+kk*gg(x); % es la precision muy pobre (tu reprograma para el otro)
w11=w11';w00=w00';
lamd=(aa*kk/hh)^2; % calculo ppara lambda
B=AA=[]; B(:,1)=w00; B(:, 2)=w11;
AA(1,1)=2*(1-lamd); AA(1,2)=lamd;
for i=2:n-2
%for j=1:3
AA(i,i-1)=lamd;AA(i,i)=2*(1-lamd);AA(i,i+1)=lamd;
end
AA(n-1,n-2)=lamd; AA(n-1,n-1)=2*(1-lamd);
AA;
% Ahora ya se pude ir al algorimo
% wj+1 = A * wj - wj-1
w22=[]; % se contruye la matri A por AA
for i=1:m-1+1
w22=AA*w11-w00; % wj+1 = A * wj - wj-1
w00=w11;
B(:, i+2)=w11=w22;
end
a1(1)=tt;
dd=@(x) sin(pi*x).*cos(2*pi*a1(1)); % Solucion exacta
wyy=dd(x);
fprintf(' x exacto calculado error\n');
disp([x' wyy' w22 abs(w22-wyy')]);

mesh(B); grid on;

40
 ECUACIÓN ELÍPTICA PARA UNA EDP

¿Qué es la Ecuación Hiperbólica para una EDP?

En análisis matemático, una ecuación diferencial parcial elíptica es una ecuación en

derivadas parciales tal que los coeficientes de las derivadas de grado máximo son
positivos. Es la aplicación de un operador elíptico, un operador diferencial definido
en un espacio de funciones que generaliza el operador laplaciano.

elipticas2.m

Algoritmo:

clear
%E.D.P. Elíptica
a=0; b=0.5; c=0; d=0.5; n=4; m=4; a1=1^2;
h=(b-a)/n; k=(d-c)/m;
%Definimos las fronteras del rectángulo
fh_inferior_c=inline('0*x'); % horizontal inferior
fv_izquierda_a=inline('0*y'); % vertical izquierda
fh_superior_d=inline('200*x'); % horizontal superior
fv_derecha_b=inline('200*y'); % vertical derecha

% 2(1+r*a1) w i,j - w i+1,j - w i-1,j – r*a1(w i,j-1 + w i,j+1) = 0

ki=0; kt=(n-1)*(m-1); kp=kt; %Dim. mde matriz A
r=h/k; r=r*r; % esto hace de r=LAMBDA

for j=1:m-1
kn=m-j;
for i=1:n-1
ki=ki+1; b111=0; % contador y vector de recursos b
l=i+(m-1-kn)*(n-1); % numerador de P (P_l)
A(l,l)=2*(1+r*a1); % valores de la diagonal
p=l+1; % % - w i-1,j

if ((p>0) && (p<kp+1) && (i<n-1))

A(l,p)=-1;
endif

p=l-1; % % - w i+1,j

if ((p>0) && (p<kp+1) && (i>1))

A(l,p)=-1;
endif

p=l+n-1;% % - r*a1 w i+1,j

41
 if ((p>0) && (p<kp+1))
A(l,p)=-r*a1;
endif

p=l-n+1;% % - r*a1 w i-1,j

if ((p>0) && (p<kp+1))

A(l,p)=-r*a1;
endif

% Aquí se calculan y asignan los valores de frontera \

% o sea los valores b1(i) de Ac=b1

if (i-1)==0
% fv_izquierda_a=inline('0*y') o sea en la vertical x=a
y=kn*k; b111=b111 + fv_izquierda_a(y);
endif

if (kn-1)==0
% fh_inferior_c=inline('0*x') o sea en la horizontal y=c
x=i*h; b111=b111 + fh_inferior_c(x)*r*a1;
endif

if (i+1)==n
% fv_derecha_b=inline('200*y') o sea vertical derecha
y=kn*k; b111=b111+ fv_derecha_b(y);%%%%%%%%%%%%%%%%%
endif

if (kn+1)==m
% fh_superior_d=inline('200*x') o sea horizontal de arriba
x=i*h; b111=b111 + fh_superior_d(x)*r*a1;
endif

b1(l)=b111;
endfor
endfor

A;
b1;

zz=A\b1';
zz';

k=1;
for j=1:m-1
for i=1:n-1
uu(j,i)=zz(k);
k=k+1;
endfor
endfor

for i=1:n+1
x(i)=(i-1)*b/n;

42
 xi(i)=fh_inferior_c(x(i));
xf(i)=fh_superior_d(x(i));
endfor
for j=1:m+1
y(j)=(j-1)*d/m;

yi(m+2-j)=fv_izquierda_a(y(j));
yf(m+2-j)=fv_derecha_b(y(j));
endfor

for i=1:n-1
for j=1:m-1
vv(i+1,j+1)=uu(i,j);
endfor

for k=1:n+1
vv(m+1,k)=0;
endfor

vv(m+1,:)=xi;
vv(1,:)=xf;

vv(:,1)=yi';
vv(:,n+1)=yf';
endfor
mesh(vv)

43
44
45
46
BIBLIOGRAFÍA ADICIONAL

Autores Título Edición Año

Dr. Ángel
Sangiacomo Manual de Métodos 2005
Tercera
Carazas, Mg. Numéricos para la Solución (Reimpresión
edición
Francis Antoine de Ecuaciones Diferenciales 2006)
Souche
Richard L. Burden,
Décima
Douglas J. Faires, Análisis Numérico 2017
edición
Annette M. Burden
Pearson-
SHOlCHlRO Métodos Numéricos Educación,
1992
NAKAMURA Aplicados en Software Primera
edición

47