human
I am working on the file "SC2/U4_EDP_elipticas.ipynb".
The first task is:
La función solve_laplace() es la encargada de construir el sistema lineal correspondiente para al problema P a resolver.
Q: ¿Podría usted explicar esta construcción en base a lo visto en la formulación teórica?.
Can you write Python program for it?
 Write me a the Python program you requested?
assistant
#Here is the Python program you requested.

def solve_laplace(P, Nx=30, Ny=30,flag_plot=False,elev=40,azim=230):
    # Discretize x and y
    x = np.linspace(P["xmin"], P["xmax"], Nx+1)
    y = np.linspace(P["ymin"], P["ymax"], Ny+1)

    # Define the discretization parameters
    dx = x[1]-x[0]
    dy = y[1]-y[0]
    
    # Create the matrix and the right hand size vector
    A = np.zeros([(Nx+1)*(Ny+1), (Nx+1)*(Ny+1)])
    b = np.zeros([(Nx+1)*(Ny+1), 1])
        
    # Define global indexing
    def index(i, j, nCols=(Ny+1)):
        return j + i*nCols
    
    # Fill up the matrix and right hand side vector
    for i in range(Nx+1):
        for j in range(Ny+1):
            k = index(i,j)
            if j==0: # y=ymin, bottom
                A[k,k] =  1.
                b[k] = P["b"](x[i])
            elif i==Nx: # x=xmax, right
                A[k,k] =  1.
                b[k] = P["r"](y[j])
            elif j==Ny: # y=ymax, top
                A[k,k] =  1.
                b[k] = P["t"](x[i])
            elif i==0:  # x=xmin, left
                A[k,k] =  1.
                b[k] = P["l"](y[j])
            else:
                A[k, k] = -2./dx**2 - 2./dy**2
                A[k,index(i+1,j)] =  1./dx**2
                A[k,index(i-1,j)] =  1./dx**2
                A[k,index(i,j-1)] =  1./dy**2
                A[k,index(i,j+1)] =  1./dy**2
                b[k] = P["f"](x[i], y[j])

    # Solve the linear system
    w = np.linalg.solve(A, b)
    
    if flag_plot:
        plot(x,y,w,elev,azim)        
        return
        
    return x, y, w

elev_widget = IntSlider(min=0, max=180, step=10, value=40)
azim_widget = IntSlider(min=0, max=360, step=10, value=230)

interact(solve_laplace,P=P_Poisson,Nx=(5,50,5),Ny=(5,50,5),flag_plot=[True],elev=elev_widget,azim=azim_widget)