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

% 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

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)