Pascal/Turbo Pascal - Gradientes Conjugados

 
Vista:

Gradientes Conjugados

Publicado por Alan89 (4 intervenciones) el 09/05/2011 07:21:54
Buenas, queria saber si me pueden pasar un procedimiento para poder hacer gradientes conjugados??????
Me mate intentando crear uno pero la verdad es qeu entiendo poco y nada como hacerlo..

Saludos
Valora esta pregunta
Me gusta: Está pregunta es útil y esta claraNo me gusta: Está pregunta no esta clara o no es útil
0
Responder
Imágen de perfil de Jimmy Aguilar

Gradientes Conjugados

Publicado por Jimmy Aguilar (1 intervención) el 09/05/2011 18:35:12
function [x] = conjgrad(A,b,x0)

r = b - A*x0;
w = -r;
z = A*w;
a = (r'*w)/(w'*z);
x = x0 + a*w;
B = 0;

for i = 1:size(A)(1);
r = r - a*z;
if( norm(r) < 1e-10 )
break;
end if
B = (r'*z)/(w'*z);
w = -r + B*w;
z = A*w;
a = (r'*w)/(w'*z);
x = x + a*w;
end

end
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Gradientes Conjugados

Publicado por ramon (2158 intervenciones) el 11/05/2011 14:11:38
{Espero sea esto lo que pedías}


program gradientes_conjugados;
uses
crt;
const
nu = 2;
type
vect = array[1..nu] of real;
vectdoble = array[1..nu * nu] of real;
var
matriz : vectdoble;
semilla , b, ax : vect;
alfa1, x1, ax1, dk, rk1, aux2, xk : vect;
axk, rk, matrizc1, zz : vect;
delta0, alfad2, aux3 : real;
alfak, delta1, tol, betak : real;
i, fil, col, j, k, kmax : integer;



procedure producto_matriz_vector(m : vectdoble; v,ans : vect; nu : integer);
var
i, j : integer;

begin
fillchar(ans,sizeof(ans),0);
for i := 1 to nu do
for j := 1 to nu do
ans[i] := ans[i] + (m[i * nu + j] * v[j]);
end;




procedure procunto_vector_vector(v1, v2 : vect; ans2 : real; nu : integer);
var
j : integer;

begin
ans2 := 0;
for j := 1 to nu do
ans2 := ans2 + v1[j] * v2[j];
end;



procedure multiplica_componentes_vectores(v1, v2, ans : vect; nu : integer);
var
j : integer;

begin
for j := 1 to nu do
ans[j] := v1[j] * v2[j];
end;


procedure entrada_valores_matriz;
begin
writeln('Entrada Valores Matriz');
gotoxy(10,3);
for i := 1 to (nu * nu) do
begin
col := i mod nu + 1;
fil := (i - i mod nu) div nu + 1;
gotoxy(10,2 + i);write('Elemento = ',fil,' ',col);
gotoxy(28,2 + i);
readln(matriz[i]);
end;
end;

procedure entrada_vector_independiente;
begin
writeln('Entrada Valores Independiente');
gotoxy(10,3);
for i := 1 to nu do
begin
col := + 1;
fil := i;
gotoxy(10,2 + i);write('Elemento = ',fil,' ',col);
gotoxy(28,2 + i);readln(b[i]);
end;
end;

procedure entrada_semilla;
begin
writeln('Entrada Valores Semilla');
gotoxy(10,3);
for i := 1 to nu do
begin
col := + 1;
fil := i;
gotoxy(10,2 + i);write('Elemento = ',fil,' ',col);
gotoxy(28,2 + i);readln(semilla[i]);
end;
end;

procedure tolerancia;
begin
writeln('Entrada Valor Tolerancia');
gotoxy(10,3);
gotoxy(10,3);write('Tolerancia = ');
gotoxy( 23,3);readln(tol);
end;

procedure interaciones;
begin
writeln('Entrada Valor Interacion max = ',(nu * nu));
gotoxy(10,3);
gotoxy(10,3);write('Interaciones = ');
gotoxy(25,3);readln(kmax);
if kmax > (nu * nu) then
kmax := (nu * nu);
end;

procedure proceso_iniciado;
begin
for i := 1 to nu do
fil := i;
matrizc1[i] := 1 / (matriz[fil * nu + fil]);
producto_matriz_vector(matriz, semilla,ax,nu);
for i := 1 to nu do
rk[i] := b[i] - ax[i];
multiplica_componentes_vectores(matrizc1,rk, zz, nu);
procunto_vector_vector(rk, zz, delta0, nu);
producto_matriz_vector(matriz, rk, alfa1, nu);
procunto_vector_vector(rk, alfa1, alfad2, nu);
if (delta0 > 0) and (alfad2 > 0) then
alfak := delta0 / alfad2
else
alfak := 0.1;
for j := 1 to nu do
xk[j] := semilla[j] + alfak * rk[j];
producto_matriz_vector(matriz, xk, ax1, nu);
for i := 1 to nu do
rk[i] := b[i] - ax[i];
delta1 := 0;
multiplica_componentes_vectores(matrizc1, rk, zz, nu);
procunto_vector_vector(rk, zz, delta1, nu);
writeln('Iteracion Activa');
k := 0;
while (delta1 > tol) or (k < kmax) do
begin
if (delta1 > 0) and (delta0 > 0) then
betak := delta1 / delta0
else
betak := 0.1;
for i := 1 to nu do
dk[i] := zz[i] + betak * dk[i];
producto_matriz_vector(matriz, dk, aux2, nu);
procunto_vector_vector(aux2, dk, aux3, nu);
if (delta1 > 0) and (aux3 > 0) then
alfak := delta1 / aux3
else
alfak := 0.1;
for i := 1 to nu do
xk[i] := xk[i] + alfak * dk[i];
for i := 1 to nu do
rk1[i] := rk[i];
producto_matriz_vector(matriz, xk, axk, nu);
for i := 1 to nu do
rk[i] := axk[i] + b[i];
multiplica_componentes_vectores( matrizc1, rk, zz, nu);
delta0 := delta1;
procunto_vector_vector(rk, zz, delta1, nu);
k := k + 1;
end;
writeln('Numero total de Iteraciones = ',k);
for j := 1 to k do
writeln('Valor ',j,' = ',xk[j]:8:2);
readln;
end;

begin
clrscr;
textcolor(15);
entrada_valores_matriz;
clrscr;
entrada_vector_independiente;
clrscr;
entrada_semilla;
clrscr;
tolerancia;
clrscr;
interaciones;
clrscr;
proceso_iniciado;
end.
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Gradientes Conjugados

Publicado por Alan89 (4 intervenciones) el 09/05/2011 20:05:49
Muchisimas gracias!
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Gradientes Conjugados

Publicado por Alan 89 (4 intervenciones) el 12/05/2011 01:02:10
Ah... tremendo, con razon no me salia programarlo
Muchas gracias :D
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Gradientes Conjugados

Publicado por Alan89 (4 intervenciones) el 09/05/2011 20:07:32
Eso es para Math Lab, cualquiera hace copy-past del wiki....
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar