Python - Transformación de un código de Matlab a Python

   
Vista:

Transformación de un código de Matlab a Python

Publicado por Ayuda a transformar un código de MatLab a Python (1 intervención) el 21/10/2017 03:03:17
Código en Matlab

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
% clear workspaces
clear
clc
 
% Inicializacion de parametros
y0 = 10;
yf = 0;
x0 = 0;
xf = 5;
c = 0.02;
dx = 0.99;
dy = 0.99;
n = round((xf-x0)/dx)+5;
m = round((y0-yf)/dy)+4;
 
s= @(x,y) x*y ;%x*y;
sx = @(x,y) y;%exp(-(x+y)^2);
sy = @(x,y) x;%exp(-(x+y)^2);
 
hf = @(x) c*x+y0;
 
 
h=zeros(m,n);
 
 
for j = 3:n-2;
    x = (j-2)*dx;
    h(1,j) = hf(x);
    h(2,j)=hf(x);
end
 
for i=1:m;
    h(i,2)=h(i,4);
    h(i,1)=h(i,5);
    h(i,n)=h(i,n-4);
    h(i,n-1)=h(i,n-3);
end
 
 
for j=1:n;
    h(m,j)=h(m-4,j);
    h(m-1,j)=h(m-3,j);
end
 
h
numit = 0;
tol = 5e-3;
nmax = 10e+4;
err = 1.0;
 
while (1)
    numit = numit+1
 
    for i=3:m-2;
        for j=3:n-2;
            oldval = h(i,j);
            h(i,j) = ((-s(i*dx,j*dy)+dx*sx(i*dx,j*dy))*h(i+2,j)...
                   + 8*(2.0*s(i*dx,j*dy)+dx*sx(i*dx,j*dy))*h(i+1,j)...
                   + 8*(2.0*s(i*dx,j*dy)-dx*sx(i*dx,j*dy))*h(i-1,j)...
                   - (s(i*dx,j*dy)-dx*sx(i*dx,j*dy))*h(i-2,j)...
                   - (s(i*dx,j*dy)+dx*sy(i*dx,j*dy))*h(i,j+2)...
                   + 8*(2.0*s(i*dx,j*dy)+dx*sy(i*dx,j*dy))*h(i,j+1)...
                   + 8*(2.0*s(i*dx,j*dy)-dx*sy(i*dx,j*dy))*h(i,j-1)...
                   -(s(i*dx,j*dy)-dx*sy(i*dx,j*dy))*h(i,j-2))/(60.0*s(i*dx,j*dy));
            err = abs((h(i,j)-oldval)/h(i,j));
 
        end
    end
 
for j=1:n;
    h(m,j)=h(m-4,j);
    h(m-1,j)=h(m-3,j);
end
 
 
    for i=1:m;
        h(i,2)=h(i,4);  %Se satisfacen condiciones del lado izquierdo
        h(i,1)=h(i,5);
        h(i,n)=h(i,n-4); %Se satosfacen condiciones del lado derecho
        h(i,n-1)=h(i,n-3);
 
    end
     h       %identificar valores previos a la iteracion final
    err     %error de cada iteracion
     if err<tol || numit>nmax, break,end
 
end
 
% impresion de resultados y grafica
h
 
[x,y] = meshgrid(1:1:n,1:1:m);   % generador de la malla, :) inverti 7 y 13
figure
surface(dx*x,dy*y,h)   % graficador
shading interp
xlabel('x')
ylabel('y')
zlabel('h(x,y)')
title('Nivel freático de Pendiente lineal---Orden4 ')
view (3)  % controla grafico 3d



ahora quiero transformar el código anterior a Python pero ya no se que más hacer

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
import numpy as n
import math as m
import matplotlib.pyplot as p
# INICIALIZACIÓN DE PARÁMETROS
 
y0 = 10.0
yf = 0.0
x0 = 0.0
xf = 5.0
c = 0.02
dx = 0.99
dy = 0.99
nx = round((xf-x0)/dx)+5  # numero de columnas ***se le sumaron dos columnas porque es una aproximación de orden 4
ny = round((y0-yf)/dy)+4
 
# DEFINICIÓN DE SIGMA Y SUS PRIMERAS DERIVADAS
 
def s(x,y): return x*y
def sx(x,y): return y
def sy(x,y): return x
 
# DEFINICIÓN DE LA CONDICIÓN DE FRONTERA
 
hf = lambda x: c*x + y0
 
# DEFINIMOS LA MATRIZ DE CALOR
 
h=n.zeros((ny,nx))
 
# RELLENAMOS VALORES EN LA FRONTERA SUPERIOR
 
for j in range(2,nx-1):
    x = (j-1)*dx
    h[0,j] = hf(x)
    h[1,j] = hf(x)
 
# SATISFACEMOS CONDICIONES DE FRONTERA LATERALES
 
for i in range(0,ny):
    h[i,1] = h[i,3]  # Se satisfacen condiciones del lado izquierdo 
    h[i,0] = h[i,4]
    h[i,nx-1] = h[i,nx-5]
    h[i,nx-2] = h[i,nx-4]
 
# SATISFACEMOS LAS CONDICIONES DE FRONTERA INFERIOR
for j in range(0,nx):
    h[ny-1,j] = h[ny-5,j];
    h[ny-2,j] = h[ny-4,j];
 
print (h)
# PARÁMETROS DE LAS ITERACIONES
 
numit = 0  # NÚMERO DE ITERACIONES 
tol = 5e-03
nmax = 10e+04
err = 1.0
errores=n.zeros((nmax))
print(errores)
 
 
while (tol<err) and (numit<nmax):
    errores[numit]=err;
    numit += 1
    #print(numit)
    # CÁLCULO DE PUNTOS INTERIORES CON EL OPERADOR DE 9 PUNTOS
    for i in range(3,ny-2):
        for j in range(3,nx-2):
            oldval = h[i,j]
            h[i,j] = ((-s(i*dx,j*dy)+dx*sx(i*dx,j*dy))*h[i+2,j] \
                   + 8*(2.0*s(i*dx,j*dy)+dx*sx(i*dx,j*dy))*h[i+1,j] \
                   + 8*(2.0*s(i*dx,j*dy)-dx*sx(i*dx,j*dy))*h[i-1,j] \
                   - (s(i*dx,j*dy)-dx*sx(i*dx,j*dy))*h[i-2,j] \
                   - (s(i*dx,j*dy)+dx*sy(i*dx,j*dy))*h[i,j+2] \
                   + 8*(2.0*s(i*dx,j*dy)+dx*sy(i*dx,j*dy))*h[i,j+1] \
                   + 8*(2.0*s(i*dx,j*dy)-dx*sy(i*dx,j*dy))*h[i,j-1] \
                   -(s(i*dx,j*dy)-dx*sy(i*dx,j*dy))*h[i,j-2]) \
                   /(60.0*s(i*dx,j*dy))
            err = abs((h[i,j]-oldval)/h[i,j]);
    #print(err)
 
    for j in range (0,nx):
        h[ny-1,j] = h[ny-5,j]
        h[ny-2,j] =h [ny-4,j]
 
        for i in range (0,ny):
            h[i,1] = h[i,3]  #Se satisfacen condiciones del lado izquierdo 
            h[i,0] = h[i,4]
            h[i,nx-1] = h[i,nx-5]   #Se satosfacen condiciones del lado derecho
            h[i,nx-2] = h[i,nx-4]
            print(h);
            print (err)
            if err<tol or numit>nmax:
               break
 
print(h)
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