Pablo Millanao

miércoles, 18 de enero de 2012

Runge Kutta 44 en Matlab

clear all;close all;clc
fun=input('Ingresa la función f(t,y) entre comillas: ');
f=inline(fun,'t','y');
y(1)=input('Ingrese la condición inicial y(0): ');
t0=input('Ingrese t0: ');
max=input('Ingrese el t final: ');
h=input('Ingrese el paso: ');
t=t0:h:max;
n=length(t);

for i=1:n-1

k1=feval(f,t(i),y(i));
k2=feval(f,t(i)+h/2,y(i)+h*k1/2);
k3=feval(f,t(i)+h/2,y(i)+h*k2/2);
k4=feval(f,t(i)+h,y(i)+h*k3);

y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);

end

plot(t,y)

Euler explícito en Matlab

clear all;close all;clc
fun=input('Ingresa la función f(t,y) entre comillas: ');
f=inline(fun,'t','y');
y(1)=input('Ingrese la condición inicial y(0): ');
t0=input('Ingrese t0: ');
max=input('Ingrese el t máximo: ');
h=input('Ingrese el paso: ');
t=t0:h:max;
n=length(t);

for i=1:n-1
y(i+1)=y(i)+h*feval(f,t(i),y(i));
end

plot(t,y)

martes, 17 de enero de 2012

Newton Raphson Multivariable

Si quieren que el usuario ingrese las funciones, jacobiana, tolerancia, etc. deben ocupar el comando input.

%newton raphson multivariable

clear all; close all; clc;

%funciones  x^2+x*y+y^2-1=0

%           x^2-y=0

f=inline('[x(1)^2+x(1)*x(2)+x(2)^2-1; x(2)-x(1)^2]');

df=inline('[2*x(1)+x(2) x(1)+2*x(2); -2*x(1) 1]');

%punto de partida para x e y

x0=[1;1];

%iteraciones máximas

maxi=10;

%tolerancia

tol=1.e-12;

corr=tol+1;

k=0;

while k<=maxi & norm(corr,inf)>tol

    fx0=feval(f,x0);

    dfx0=feval(df,x0);

    corr=dfx0\fx0;

    x=x0-corr;

    x0=x;

    k=k+1;

end

disp('la solución es')

disp(x)

if norm(corr,inf)>tol

    error('Se alcanzó el número máximo de iteraciones')

end

disp(['se calculo en ' num2str(k) ' iteraciones'])

Bisección en Matlab

%bisección
clear all; close all; clc;
fun=input('Ingrese la función f(x) entre comillas: ');
f=inline(fun);
a=input('Ingrese el límite inferior del intervalo: ');
b=input('Ingrese el límite superior del intervalo: ');
tol=input('Ingrese la tolerancia: ');
maxit=input('Ingrese el máximo de iteraciones: ');

fa=feval(f,a);
fb=feval(f,b);
k=0;

if fa*fb>0
error('la función tiene el mismo signo en ambos extremos')
end

while abs(b-a)>tol & k<maxit
if fa*fb==0

if fa==0
s=a;
break
else
s=b;
break
end

else

k=k+1;
s=(a+b)/2;
fs=feval(f,s);
if fa*fs<0
b=s;
fb=feval(f,b);
else
a=s;
fb=feval(f,a);
end
end
end

disp('la solución es')
disp(s)

if abs(b-a)>tol
error('Se alcanzó el número máximo de iteraciones')
end

disp(['se calculo en ' num2str(k) ' iteraciones'])

lunes, 12 de diciembre de 2011

Ecuaciones diferenciales ordinarias

'edo a resolver
Public Function g(t, y) As Double

g = (t - y) / 2

End Function

'runge kutta 44
Public Function rk_44(t0, y0, h) As Double

k1 = g(t0, y0)
k2 = g(t0 + h / 2, y0 + h * k1 / 2)
k3 = g(t0 + h / 2, y0 + h * k2 / 2)
k4 = g(t0 + h, y0 + h * k3)

rk_44 = y0 + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)

End Function

'heun
Public Function e_heun(t0, y0, h) As Double

yp = y0 + h * g(t0, y0)
e_heun = y0 + 0.5 * h * (g(t0, y0) + g(t0 + h, yp))

End Function

'poligono mejorado
Public Function poli_m(t0, y0, h) As Double

ym = y0 + 0.5 * h * g(t0, y0)
poli_m = y0 + h * g(t0 + h / 2, ym)

End Function

'ralston
Public Function ralston_(t0, y0, h) As Double

k1 = g(t0, y0)
k2 = g(t0 + 0.75 * h, y0 + 0.75 * h * k1)
ralston_ = y0 + (h / 3) * (k1 + 2 * k2)

End Function

'euler explicito
Public Function eu_ex(t0, y0, h)

eu_ex = y0 + h * g(t0, y0)

End Function

Invertir Matrices a través de VBA

Podemos calcular una matriz de 2x2 y de 3x3 en VBA con una subrutina a partir de la regla

$$A^{-1}=\dfrac {1} {det\left( A\right) }adj\left( A\right)$$

a continuación se muestran los códigos, para que esta subrutina corra de forma correcta la matriz a invertir debe encontrarse emplazada en el rango de celdas "B3:D5" y los resultados serán desplegados en el rango de celdas "F3:H5" como se muestra en la figura

Private Sub inv3x3()

'invertir una matriz A = [a b;c d]
'recordando que la la inversa de esta matriz es:
'1/det(A).*[d -b;-c a]

Dim A(1 To 2, 1 To 2) As Double

For i = 1 To 2
For j = 1 To 2

A(i, j) = Cells(2 + i, 1 + j)

Next j
Next i

'determinante de la matriz
det = A(1, 1) * A(2, 2) - A(1, 2) * A(2, 1)

'si el determinate de la matriz es 0 => inversa de la matriz no existe
If det = 0 Then
MsgBox ("La matriz no es invertible")
Exit Sub
End If

'reordenando los terminos y multiplicando por el reciproco del determinate
Cells(3, 5) = det ^ -1 * A(2, 2)
Cells(3, 6) = det ^ -1 * -A(1, 2)
Cells(4, 5) = det ^ -1 * -A(2, 1)
Cells(4, 6) = det ^ -1 * A(1, 1)

End Sub

para invertir una matriz de 2x2 se ocupa el siguiente código, para que esta subrutina se ejecute de forma optima la matriz a invertir debe ubicarse en el rango "B3:C4" y los resultados serán arrojados en rango de celdas "E3:F4" como se muestra en la siguiente imagen

Private Sub inv2x2()

'invertir una matriz A = [a b;c d]
'recordando que la la inversa de esta matriz es:
'1/det(A).*[d -b;-c a]

Dim A(1 To 2, 1 To 2) As Double

For i = 1 To 2
For j = 1 To 2

A(i, j) = Cells(2 + i, 1 + j)

Next j
Next i

'determinante de la matriz
det = A(1, 1) * A(2, 2) - A(1, 2) * A(2, 1)

'si el determinate de la matriz es 0 => inversa de la matriz no existe
If det = 0 Then
MsgBox ("La matriz no es invertible")
Exit Sub
End If

'reordenando los terminos y multiplicando por el reciproco del determinate
Cells(3, 5) = det ^ -1 * A(2, 2)
Cells(3, 6) = det ^ -1 * -A(1, 2)
Cells(4, 5) = det ^ -1 * -A(2, 1)
Cells(4, 6) = det ^ -1 * A(1, 1)

End Sub

domingo, 11 de diciembre de 2011

Ecuaciones no lineales

Funciones programadas en VBA para resolver ecuaciones no lineales

'función a resolver
Public Function g(x) As Double
g = Exp(-x) - x
End Function

'derivada de g(x)
Public Function dg(x) As Double
h = 0.000001
dg = (g(x + h) - g(x)) / h
End Function

'bisección

Public Function biseccion(a, b, n)

If g(a) * g(b) > 0 Then
MsgBox ("la función tiene el mismo signo en ambos extremos")
End If

If g(a) * g(b) = 0 Then

If g(a) = 0 Then
s = a
Else
s = b
End If

Else

For i = 1 To n

s = (a + b) / 2

If g(a) * g(s) < 0 Then
b = s
Else
a = s
End If

Next i

End If

biseccion = s

End Function

'newton raphson
Public Function newtonraphson(x0, n) As Double

xi=x0

For i = 1 To n

xi_1 = xi - g(xi) / dg(xi)
xi = xi_1

Next i

newtonraphson = xi_1

End Function

'punto fijo
Public Function ptofijo(x0, n)

xi=x0

For i = 1 To n

x_i = g(xi) + xi
xi = xi_1

Next i

ptofijo = x_i

End Function