Algoritmos de Optimización No Lineal

Descargar como pdf o txt
Descargar como pdf o txt
Está en la página 1de 31

INFORME FINAL

OPTIMIZACIÓN NO LINEAL

PROFESOR: DOCTOR ANGEL SANGIACOMO CARAZAS

ALUMNO: LARICO ALVIZ JOEL NELSON

19 DE AGOSTO DE 2022
INFORME SOBRE LOS ALGORITMOS DE OPTIMIZACIÓN NO LINEAL

REDACTADO POR

Larico Alviz Joel Nelson

PROGRAMAS Y/O HERRAMIENTAS UTILIZADAS

 Octave
 Calculadora Científica
 Libro de Optimización no Lineal

OBJETIVOS GENERALES

 Encontrar soluciones (a través de métodos matemáticos con el uso de sistemas lineales)


a problemas de carácter económico-técnico, representados por la limitación de
recursos.
 Resolver casos de combinación óptima de mezclas de producción, disposición interna
de procesos, maximización de beneficios, localización, asignación de recursos,
minimización de costos, transporte, entre otros.

ALGORITMOS DE PROGRAMACIÓN

La programación no lineal es un método por el cual se optimiza, ya sea maximizando


o minimizando, una función objetivo. Esto, tomando en cuenta distintas restricciones
dadas. Se caracteriza porque la función objetivo, o alguna de las restricciones, pueden
ser no lineales calculados por diferentes programas mediante distintos algoritmos.
ALGORITMOS APLICADOS EN CLASE

MÉTODO DE LAS CUERDAS

¿Qué es el Método de las Cuerdas?

En cálculo numérico, el método de la regula falsi (regla del falso) o falsa posición es
un método iterativo de resolución numérica de ecuaciones no lineales. El método
combina el método de bisección y el método de la secante.

Como en el método de bisección, se parte de un intervalo inicial [a0,b0] con f(a0) y


f(b0) de signos opuestos, lo que garantiza que en su interior hay al menos una raíz
(véase Teorema de Bolzano). El algoritmo va obteniendo sucesivamente en cada paso
un intervalo más pequeño [ak, bk] que sigue incluyendo una raíz de la función f.

A partir de un intervalo [ak, bk] se calcula un punto interior ck:


3
Dicho punto es la intersección de la recta que pasa por (a,f(ak)) y (b,f(bk)) con el eje
de abscisas (igual a como se hace en el método de la secante).

Se evalúa entonces f(ck). Si es suficientemente pequeño, ck es la raíz buscada. Si no,


el próximo intervalo [ak+1, bk+1] será:

 [ak, ck] si f(ak) y f(ck) tienen signos opuestos;


 [ck, bk] en caso contrario.

Se puede demostrar que bajo ciertas condiciones el método de la falsa posición tiene
orden de convergencia lineal, por lo que suele converger más lentamente a la
solución de la ecuación que el método de la secante, aunque a diferencia de en el
método de la secante el método de la falsa posición siempre converge a una solución
de la ecuación.

 cuerdas_11.m

Algoritmo:

f=@(x) 4.5 - x.^2;

a=2; b=3; t0=a; ert=0.001;


for i=1:10
t= a - f(a)*(b-a)/(f(b)-f(a));
fprintf('%d a=%f t=%f b=%f',i, a, t, b);
if f(t)*f(a)<0 b=t;
else a=t;
end;
a1=abs(t-t0); t0=t;
fprintf(' er=%f\n',a1);
if a1<ert break; end;

end;

%fprintf('%d a=%f t=%f b=%f Er=%f\n',i, a, b, )

MÉTODO DE LOS MÍNIMOS LOCALES

 ¿Qué es el Método de los Mínimos Locales?

4
En matemáticas, los máximos y mínimos de una función, conocidos colectivamente
como extremos de una función, son los valores más grandes (máximos) o más
pequeños (mínimos), que toma una función en un punto situado ya sea dentro de
una región en particular de la curva (extremo local o relativo) o en el dominio de la
función en su totalidad (extremo global o absoluto).123 De manera más general, los
máximos y mínimos de un conjunto (como se define en teoría de conjuntos) son
los elementos mayor y menor en el conjunto, cuando existen. El localizar valores
extremos es el objetivo básico de la optimización matemática.

 minimo_11.m

Algoritmo:

%% EL DE LA CLSE CON CLA FLAQUITA?????????

ff=@(x,y) 3*x.^2-4*x.*y+x-3*y+4*y.^2+3;

% algritmo
clc;
x0=[3, 8];
for j=1:5
d=0;
xx=x0-[d,0];
k1=ff(xx(1),xx(2))
for i=1:5
d=d-.05;
xx=x0-[d,0]
k2=ff(xx(1),xx(2))
if k1<k2 break; end;
k1=k2;
x0=xx;
end;
5
xx
disp('===========================');
%%%%%%%%%%%%% (0,1)
d=0;
xx=x0-[0,d];
k1=ff(xx(1),xx(2))
for i=1:5
d=d+.05;
xx=x0-[0,d]
k2=ff(xx(1),xx(2))
if k1<k2 break; end;
k1=k2;
x0=xx;
end;
xx
disp('RRRRRRRRRRRRRRRRRRRRRRRRRRR');
end;

 minimo_22.m

Algoritmo:

% ES CORREGIDO YA CON direcciones coordenadas


% (1, 0) <==> como gradiente de x
% (0, 1) <==> como gradiente de y

ff=@(x,y) 3*x.^2-4*x.*y+x-3*y+4*y.^2+3;

AA=[6 -4; -4 8]; bb=-[1 -3]';


cc=AA\bb
% algritmo
clc;
x0=[3, 8];
for j=1:1500
d=0;
xx=x0-[d,0];yy=xx;
k1=ff(xx(1),xx(2));
for i=1:50
d=d+.005;
xx=x0-[d,0];
k2=ff(xx(1),xx(2));
if k1<k2 break; end;
k1=k2;
x0=xx;
end;
xx;
%disp('===========================');
%%%%%%%%%%%%% (0,1)
d=0;
xx=x0-[0,d];
k1=ff(xx(1),xx(2));
6
for i=1:50
d=d+.005;
xx=x0-[0,d];
k2=ff(xx(1),xx(2));
if k1<k2 break; end;
k1=k2;
x0=xx;
end;
%xx
%disp('RRRRRRRRRRRRRRRRRRRRRRRRRRR');
if norm(xx-yy) <0.001 break; end;
yy=xx;

end;
yy
k1
j

 minimo_aa11.m

Algoritmo:

% ES CORREGIDO YA CON TODO Y GRADIENTE EXACTO


% 6x - 4y + 1 <==> gradiente de x
% -4x + 8y - 3 <==> gradiente de y

ff=@(x,y) 3*x.^2-4*x.*y+x-3*y+4*y.^2+3;

AA=[6 -4; -4 8]; bb=-[1 -3]';


cc=AA\bb % solucion exacta
% algritmo
clc;
x0=[3, 8];
for j=1:1500
d=0;
xx=x0-[d,0];yy=xx;
k1=ff(xx(1),xx(2));
for i=1:5
d=d+.05;
xx=x0-d*[6*x0(1)-4*x0(2)+1,0];
k2=ff(xx(1),xx(2));
if k1<k2 break; end;
k1=k2;
x0=xx;
end;
xx;
%disp('===========================');
%%%%%%%%%%%%% (0,1)
d=0;
xx=x0-[0,d];
k1=ff(xx(1),xx(2));
7
for i=1:5
d=d+.05;
xx=x0-d*[0,-4*x0(1)+8*x0(2)-3];
k2=ff(xx(1),xx(2));
if k1<k2 break; end;
k1=k2;
x0=xx;
end;
%xx
%disp('RRRRRRRRRRRRRRRRRRRRRRRRRRR');
if norm(xx-yy) <0.00001 break; end;
yy=xx;

end;
yy
k1
j

INTERPOLACIONES CÚBICAS

 ¿Qué son las Interpolaciones Cúbicas?

En este caso, cada polinomio 𝑃(𝑥) a través del que construimos los Splines en [𝑚, 𝑛]
tiene grado 3. Esto quiere decir, que va a tener la forma 𝑃(𝑥) = 𝑎𝑥 3 + 𝑏𝑥 2 + 𝑐𝑥 + 𝑑
En este caso vamos a tener cuatro incógnitas por cada intervalo (𝑎, 𝑏, 𝑐, 𝑑), y una
nueva condición para cada punto común a dos intervalos, respecto a la derivada
segunda:
 Que las partes de la función a trozos 𝑃(𝑥) pasen por ese punto. Es decir, que las
dos 𝑃𝑛 (𝑥) que rodean al 𝑓(𝑥) que queremos aproximar, sean igual a 𝑓(𝑥) en cada
uno de estos puntos.
 Que la derivada en un punto siempre coincida para ambos "lados" de la función
definida a trozos que pasa por tal punto común.
 Que la derivada segunda en un punto siempre coincida para ambos "lados" de la
función definida a trozos que pasa por tal punto común.
Como puede deducirse al compararlo con el caso de splines cuadráticos, ahora no
nos va a faltar una sino dos ecuaciones (condiciones) para el número de incógnitas
que tenemos.
La forma de solucionar esto, determina el carácter de los splines cúbicos. Así,
podemos usar:
 Splines cúbicos naturales: La forma más típica. La derivada segunda de 𝑃 se hace
0 para el primer y último punto sobre el que está definido el conjunto de Splines,
esto son, los puntos 𝑚 y 𝑛 en el intervalo [𝑚, 𝑛].
 Dar los valores de la derivada segunda de 𝑚 y 𝑛 de forma "manual", en el
conjunto de splines definidos en el intervalo [𝑚, 𝑛].
 Hacer iguales los valores de la derivada segunda de 𝑚 y 𝑛 en el conjunto de
splines definidos en el intervalo [𝑚, 𝑛].
8
 Splines cúbicos sujetos: La derivada primera de 𝑃 debe tener el mismo valor que
la derivada primera de la función para el primer y último punto sobre el que está
definido el conjunto de Splines, esto son, los puntos 𝑚 y 𝑛 en el intervalo [𝑚, 𝑛].

 cubicas_11_a1.m

Algoritmo:

% cubicasa

f=@(x) x.*x - 5*x +1;


df=@(x) 2*x-5;

x=-3; d=1;
al=0;
k1=f(x+al*d), k2= df(x+al*d);
al1=1;
for i=1:10
al=al+1;
k3=df(x+al*d);
if k2*k3<0 break; end;

end;
a=x+(al-1)*d; b=x+al*d;

z=(3*(f(a)-f(b)))/(b-a)+(df(a)+df(b))
w=sqrt(z.*z - df(a)*df(b))
ac=b-(df(b)-z-w)/(df(b-df(a)+2*w))*(b-a)

MÉTODO DE LA GRADIENTE

 ¿Qué es el Método de la Gradiente?

9
Para la optimización de modelos de Programación No Lineal sin restricciones se
dispone de una categoría de métodos denominados «Algoritmos Generales de
Descenso» entre los cuales se destaca el Método del Gradiente o Método del
Descenso más Pronunciado (conocido adicionalmente como Método de
Cauchy) que reducen el cálculo de un mínimo local a una secuencia de problemas de
búsqueda lineal (o búsqueda unidimensional).
Consideremos el siguiente problema de Programación No Lineal no restringido:

Es importante observar lo siguiente: El punto de partida para comenzar las


iteraciones es arbitrario y al ser evaluado en la función objetivo se alcanza un valor
de V(8,7)=-149.
Si evaluamos la coordenada que se alcanza al realizar una iteración del método la
función objetivo obtiene el siguiente valor V(12,5)=-169 que como se puede apreciar
reduce el valor de la función objetivo.
En resumen el Método del Gradiente consta de 2 pasos principales:
10
Primero: El cálculo de una dirección de descenso que esta dado por el negativo del
gradiente de la función objetivo evaluado en el punto de partida o en el de la k-
ésima iteración. En el ejemplo dicha dirección desde la coordenada
original x°=(8,7) esta dada en la dirección del vector d°=(8,-4).
Segundo: Obtener la magnitud del paso α (alfa) que determina cuánto se avanza en
una determinada dirección de descenso. Esto se logra generando una función
unidimensional en términos de este parámetro (respecto a la función objetivo
original). En el ejemplo dicha magnitud del paso es α=1/2.
Finalmente se actualiza la coordenada según lo descrito previamente
alcanzando (x1,x2)=(12,5) que como se corroboró otorga un valor en la función
objetivo menor al punto de partida (arbitrario).
¿Cómo determinar si se ha alcanzado la solución óptima del problema no
restringido a través del Método del Gradiente?
Para ello se debe verificar que la dirección de descenso en la k-ésima iteración es nula
(cero).
Se puede corroborar en este ejemplo que esto se logra al intentar realizar una nueva
iteración a partir de (x1,x2)=(12,5).
Adicionalmente se generan las condiciones de segundo orden calculando la Matriz
Hessiana o de segundas derivadas de la función objetivo:

Donde los determinantes son estrictamente mayores a


cero: D1=2>0 y D2=4>0 (la Matriz Hessiana por tanto es Positiva Definida) por lo
cual la función objetivo es estrictamente convexa y dada la condición anterior es
suficiente para garantizar que (x1,x2)=(12,5) es la solución óptima del problema.

 gradiente_11.m

Algoritmo:

% gradiente
% f(x,y) = 2x^2-3xy+7y^2-3x+4y+4

ff=@(x,y) 2*x.^2-3*x.*y+7*y.^2-3*x+4*y+4;
ffx=@(x,y) 4*x-3*y-3;
ffy=@(x,y) -3*x+14*y+4;
x0=1; y0=2; er11=.001;

for i=1:20
alfa=0; k1=ff(x0,y0);
for j=1:10
alfa=alfa+.1;
11
xtx=x0-alfa*ffx(x0,y0);
xty=y0-alfa*ffy(x0,y0);
k2=ff(xtx,xty);
if k1>k2 break; end;
k1=k2;
end; % fin j
k3=norm([x0 y0]-[xtx xty]);
printf('%d x=%f y=%f Er=%f k=%f\n', i, x0, y0, k3, ff(xtx, xty));
if k3<er11 break; end;
x0=xtx; y0=xty;
end;
printf('%d x=%f y=%f Er=%f\n', i+1, xtx, xty, k3);

MÉTODO DE LA GRADIENTE

 ¿Qué es el Método de la Gradiente?


Si recordamos la historia en busca del concepto de divina proporción. Leonardo
Pisano, también conocido como Fibonacci, fue un famoso matemático de Italia que
se dedicó a divulgar por Europa el sistema de numeración árabe (1, 2, 3…) con base
decimal y con un valor nulo (el cero) en su Libro del ábaco en 1202.
Pero, el gran descubrimiento de este matemático fue la Sucesión de Fibonacci que,
posteriormente, dio lugar a la proporción áurea en arte.
Cálcular la proporcionalidad áurea:
Una herramienta de utilidad para obtener las medidas de forma rápida y práctica es
la siguiente calculadora de proporciones áureas que nos ayudará a encontrar las
medidas:
Golden Rectangle Calculator:

12
 aurea_11_a1.m

Algoritmo:

f = @(x) 2*x*x-3*x-4;
a = -3; b = 4; error = 0.001; Tau = (sqrt(5)-1)/2;
if a > b x1 = a; a = b; b = x1 end;
x2 = a + (b-a) * Tau;
x1 = a + (b-a) * (1 - Tau);
a1 = a; b1 = b; i1=1;
fx1 = f(x1); fx2 = f(x2);
do
if fx1 < fx2
b1 = x2; x2=x1; fx2=fx1;
x1 = a1 + (b1-a1) * (1-Tau); fx1 = f(x1);
else %{if fx1 > fx2 then}
a1 = x1; x1=x2; fx1=fx2;
x2 = a1 + (b1-a1) * (Tau);
fx2 = f(x2);
end;
13
i1 = i1 + 1;
until abs(a1-b1) < error;
%Str(a1:1:10,sx1);Str(b1:1:10,sx2);
e1=abs(a1-b1);
fprintf(' El intervalo Final es [ %f ; %f ] Er=%f\n', a1, b1, e1);
%str( error * (b-a):1:8 , sx1); str( b1-a1:1:8 , sx2);

%end;

 Aurea01_20220512_a11.m
Algoritmo:

%procedure TFibonacci.Button1Click(Sender: TObject);


%var a, b, x1, x2, x3, x4, Tau, error, a1, b1 : Double;
% fx1, fx2 : double;
% s, i1 : integer;
% sx1, sx2, sx3, sx4, sx5, sx6 : string[30];
f = @(x) 2*x.*x-3*x-4;
a = -2;
b = 5;
error1 = 0.001;
Tau := (sqrt(5)-1)/2;
if a > b x1 = a; a = b; b = x1 end;
x2 = a + (b-a) * Tau;
x1 = a + (b-a) * (1 - Tau);
a1 = a; b1 = b; i1 = 0; s = 3;
do %repeat
if i1 > 0

if (s=0) x2 = a1 + (b1-a1) * Tau;


if (s=1) x1 = a1 + (b1-a1) * (1 - Tau);
end; % fin if
fx1 = f(x1); fx2 = f(x2);
%str(a1:1:8, sx1); str(x1:1:8, sx2); str(x1:1:8, sx5); str(b1:1:8,
sx6);
%str(fx1:1:8, sx3); str(fx2:1:8, sx4);
%ListBox1.Items.Add(' a= '+sx1+' b= '+sx2+ ' x1= '+ sx5 + '
x2=' + sx6 +
% ' f(a)= '+sx3+ ' f(b)= '+sx4);
%(* if x1 > x2 then
% begin x4 := x1; x1 := x2; x2 := x4 end;*)
% fx1 := f(x1); fx2 := f(x2);
if fx1 > fx2
i1 = i1 + 1;
if x1 > x2
b1 = x1; s = 1;
else
a1 = x1; s = 0; x1 = x2;
end;
else %{if fx1 < fx2 then}
i1 = i1 + 1;
if x1 > x2
14
a1 = x2; s = 0;
else
b1 = x2; s = 1; x2 = x1;
end;
%(* if fx1 = fx2 then
% begin
% if x1 <> x2 then
% begin
% b1 := x2; a1 := x1; s := 5; i1 := i1 + 1; x1:= x2
% end
% else
% begin s :=1 ;
% end
% end; *)
Until abs(a1-b1) < error1;
printf(' El intervalo Final es [ %f ; %f ]\n', a1, x1);

end;

%end.

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…

 Fibonacci_11_a1.m

Algoritmo:
15
%function f( x : double ) : _r;
%begin f := x*x-3*x end;
clear all;
f =@(x) x*x-3*x;
n=0; error1=0.0001;
%Procedure Fibonacci( var n : _i); var xx : _r;
%begin Fibonaccin(1) := 1; Fibonaccin(2) := 1; n := 2;
% repeat n := n + 1;
% Fibonaccin(n) := Fibonaccin(n-1)+Fibonaccin(n-2);
% until Fibonaccin(n) > 1/error;
%end;
Fibonaccin=[];
Fibonaccin(1) = 1; Fibonaccin(2) = 1; n = 2;
do n = n + 1;
Fibonaccin(n) = Fibonaccin(n-1)+Fibonaccin(n-2);
until Fibonaccin(n) > 1/error1;
nn=n;
a = -4;
b = 5;
%error1 := 0.01;
%Fibonacci( nn );
%Memo1.Clear;
if a > b x1 = a; a = b; b = x1; end;
hh = (b-a) * Fibonaccin(nn-2) / Fibonaccin(nn);
x1=a+hh; x2=b-hh; fx1 = f(hh+a); fx2 = f(-hh+b); i1=nn;nn=nn;
do %nn=nn-1;
%sx3:=''; str(a:1:5,sx1);sx3:=sx3+' a= '+sx1;str(b:1:5,sx1);sx3:=sx3+' b=
'+sx1;
% str(x1:1:5,sx1);sx3:=sx3+' x1= '+sx1;str(x2:1:5,sx1);sx3:=sx3+' x2=
'+sx1;
% str(fx1:1:5,sx1);sx3:=sx3+' fx1= '+sx1;str(fx2:1:5,sx1);sx3:=sx3+'
fx2= '+sx1;
% sx3:=sx3+' '+sx2;
% Memo1.Lines.Add(sx3);
if fx1 < fx2
if x1 < x2 %then
b = x2; %end;
else a = x2; end;
hh = (b-a) * Fibonaccin(nn-2) / Fibonaccin(nn);
x2 = b-hh; fx2 = f(x2); %sx2:=' <<';
else
%begin
if x1 < x2 a = x1;
else b = x1; end;
x1 =x2;fx1=fx2;
hh = (b-a) * Fibonaccin(nn-2) / Fibonaccin(nn);
x2 = b - hh; fx2 = f(x2); % sx2:=' >>';
end; %end;
nn=nn-1;
printf(' El intervalo Final es [ %f ; %f ] \n',a ,b );
until 3 > nn;
% sx3:='';str(a:1:6,sx1);sx3:=sx3+' ( '+sx1+' , ';

16
% str(b:1:6,sx1);sx3 := sx3 + sx1 + ' )';
% Memo1.Lines.Add(' El intervalo Final es '+sx3+ ' n= '+IntToStr(i1));
printf(' El intervalo Final es [ %f ; %f ] n=%d \n',a ,b , i1);
% str( error * (b-a):1:8 , sx1); str( b-a:1:8 , sx2);
% Memo1.Lines.Add(' Con Error Pedido de Err = ' + sx1 +
% ' y Error de intervalo = '+ sx2);
%end;

%end.

MÉTODO DE MONTECARLO

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

El método de Montecarlo1 es un método no determinista o estadístico numérico,


usado para aproximar expresiones matemáticas complejas y costosas de evaluar con
exactitud. El método se llamó así en referencia al Casino de Montecarlo (Mónaco) por
ser “la capital del juego de azar”, al ser la ruleta un generador simple de números
aleatorios. El nombre y el desarrollo sistemático de los métodos de Montecarlo datan
aproximadamente de 1944 y se mejoraron enormemente con el desarrollo de
la computadora u ordenador.
El uso de los métodos de Montecarlo como herramienta de investigación proviene
del trabajo realizado en el desarrollo de la bomba atómica durante la Segunda Guerra
Mundial en el Laboratorio Nacional de Los Álamos en EE. UU. Este trabajo
conllevaba la simulación de problemas probabilísticos
de hidrodinámica concernientes a la difusión de neutrones en el material de fisión.
Esta difusión posee un comportamiento eminentemente aleatorio. En la actualidad es
parte fundamental de los algoritmos de raytracing para la generación de imágenes
3D.

17
En la primera etapa de estas investigaciones, John von Neumann y Stanislaw
Ulam refinaron esta ruleta y los métodos "de división" de tareas. Sin embargo, el
desarrollo sistemático de estas ideas tuvo que esperar al trabajo de Harris y Herman
Kahn en 1948. Aproximadamente en el mismo año, Enrico Fermi, Nicholas
Metropolis y Ulam obtuvieron estimadores para los valores característicos de
la ecuación de Schrödinger para la captura de neutrones a nivel nuclear usando este
método.
El método de Montecarlo proporciona soluciones aproximadas a una gran variedad
de problemas matemáticos posibilitando la realización de experimentos con
muestreos de números pseudoaleatorios en una computadora. El método es aplicable
a cualquier tipo de problema, ya sea estocástico o determinista. A diferencia de
los métodos numéricos que se basan en evaluaciones en N puntos en un espacio M-
dimensional para producir una solución aproximada, el método de Montecarlo tiene
1
un error absoluto de la estimación que decrece como en virtud del teorema del
√𝑁
límite central.

 monte_carlo11.m

Algoritmo:

function y=monte_carlo11
xk=2;
for j=1:100
h=4.5 - xk*xk;
if h<0 i=-1; else i=1; end;
18
h=abs(h);
kk=0; % contador
if h<1 h=h/2;
else do kk=kk+1; until kk*kk>h;
h=kk-1; end;
aleat=rand^2;
vx=xk+i*h*aleat;
%fprintf('carlo=%f er=%f\n', vx, abs(4.5-vx*vx));
if vx<0 vx=abs(vx); end;
if abs(vx-xk) <.0000001 break; end;
xk=vx;

end;
fprintf('%d carlo=%f er=%f\n',j, vx, abs(4.5-vx*vx));

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:

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 .

19
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;
if p4<0.00001 break; end;
end;
ys
plot(ys);
grid;
end; %fin de disparo ...

function y2=f1(t, y, w)
y2=w;
20
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;

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.

21
 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));
RR(n-1)=-h*h*rr(a+h*(n-1))+(1-h*pp(a+h*(n-1))/2)*yb;
if k<1
AA
RR
22
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
end
format long

23
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
AA
RR
24
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;

MÉTODO DE LAS BASES

 base_111.m

Algoritmo:

clc
%A=[1 0 0 0;1 1+2/3 1/3*3+2*1/3 1/3+3*1/9; 1 1+2/3 2*2/3+2*2/3 2*2*2/3+3*4/9; 1
2 3 4]
A=[1 0 0 0
1 1+1/3 1/9+2/3 1/27+3/9
1 1+2/3 4/9+4/3 8/27+12/9
1 2 3 4 ]

b=[0 1/3 2/3 1]'


format rat; %escribe la soluycion en fraccion
c=A\b
x=0:.01:1;
ys=c(1)+c(2)*x+c(3)*x.^2+c(4)*x.^3;
yt=exp(-x)+x-1;
plot(x, ys, x, yt);

 base_222.m

25
Algoritmo:

%este no sale pipipipipipipipiiii


A=[1 0 0 0
1 1+1/3 1/9+2/3+2 1/27+3/9+6/3
1 1+2/3 4/9+4/3+2 8/27+12/9+12/3
0 1 0 0]

b=[0 1/3 2/3 1]'

format rat;
c=A\b
format;
x=0:.01:1;
ys=c(1)+c(2)*x+c(3)*x.^2+c(4)*x.^3;
yt=sinh(2*x)/4*sinh(2) - x/4;
yp=5*sinh(2*x)/8 - x/4;
plot(x, ys, x, yt, x, yp);
grid;

 base_primer11.m

Algoritmo:

%debe salir graficaaaaaa

function yta=base_primer11
a=0; b=1; n=3; h=(b-a)/n;

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

%Construimos b está arriba


format rat;
c=A\b'
format;
x=0:.1:1;

26
ys=c(1)+c(2)*x+c(3)*x.^2+c(4)*x.^3;

yt=exp(-x)+x-1;

plot(x, ys, x, yt);


grid;

%x_x(2,3)

% 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),
end;

function v=rr(t)
v=t;

 base_primer12.m

Algoritmo:

%PROGRAMA PARA
%edo
function yta=base_primer12
a=0; b=1; n=6; h=(b-a)/n;

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

%Construimos b está arriba


format rat;

27
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(c)
k=0;
for j=1:length(c)
k=k+c(j)*x(i).^(j-1);
endfor
ys(i)=k;
end;

yt=exp(-x)+x-1; yt=[];
plot(x, ys, x, yt);
grid;

%x_x(2,3)

% 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),
end;

function v=rr(t)
v=t;

 base_primer22.m

Algoritmo:

%PROGRAMA PARA
%edo
function yta=base_primer22
a=0; b=1; n=6; h=(b-a)/n;

x1=a:h:b;
A=[]; b=[];
%contruccion 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));

28
end;
for i=2:n+1
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

%Construimos b está 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(c)
k=0;
for j=1:length(c)
k=k+c(j)*x.^(j-1);
endfor
ys(i)=k;
end;

yt=exp(-x)+x-1;
plot(x, ys, x, yt);
grid;

%x_x(2,3)

% 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),
end;

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;

29
OTROS PROGRAMAS

 comandos.m

Algoritmo:

format rat
A=[62/9 -25/9; -25/9 62/9]
b=[-1/9 -1/18]'

A\b %-0.023

disp('======================')

(sinh(2/3)/(4*sinh(2)))-(1/3)/4 %-0.033

 programa20220623.m

Algoritmo:

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(Ip1,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(1, i+1)=-IpN+Ib; A(i+1, i)=A(i, i+1);


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; h1=h;
Ic=@(z) (z-x(i-1)).*r(z); h1=x(i-1)-x(i-2); Ic=quad(Ic,x(i-2),x(i-1))/h1;
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;

30
A
b
%c = A\b'

BIBLIOGRAFÍA ADICIONAL

Autores Título Edición Año

Dr. Ángel INTRODUCCIÓN A LA


Sangiacomo OPTIMIZACIÓN NO
Arequipa
Carazas, Ing. LINEAL CON SOFTWARE Cuarta edición
2006
Fernando Ayres DE APLICACIÓN
Hidalgo DIDÁCTICA
David G. Luenberger, Linear and Nolinear Third Edition
2008
Yinyu Ye Programming (Springer)

31

También podría gustarte