<jupyter_start><jupyter_text>Singular Pertubations Problem<jupyter_code>#%matplotlib notebook
%matplotlib inline

# environment setting, before any codes
import numpy as np
import scipy.linalg as slinalg
import matplotlib.pyplot as plt
import numpy as np
import numpy.polynomial.legendre as npleg
from mpl_toolkits import mplot3d
from mpl_toolkits.mplot3d import axes3d
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
from IPython.display import clear_output, display
#The following functions will give us the exact solution of our boundary value problem

def exact_fun(x, alpha, beta, epsilon):
    return alpha + x + (beta - alpha - 1)*((np.exp(x/epsilon)-1)/(np.exp(1/epsilon)-1))

alpha = 1
beta = 3
epsilon = 0.01

number = 59
def draw_fun(number, alpha, beta, epsilon):
    x = np.linspace(0, 1, number)
    y = exact_fun(x, alpha, beta, epsilon)
    plt.plot(x,y,'r',linewidth=2)
    return epsilon
draw_fun(number, alpha, beta, epsilon)

#As we have a singular value problem, we must design a fine grid that will fit with our case

alpha = 1
beta = 3
epsilon = 0.01
m = 100
#We generate the uniform grid on [alpha, beta]

def generate_grid_y(left, right, m):
    h = (right - left)/(m+1)
    y = np.zeros(m+2)
    for j in range(m+2):
        y[j] = j*h+1
    return  y

y = generate_grid_y(alpha,beta, m)
#print(y)

#The method we use to design the grid is the dichotomy

def antecedent_k(exact_fun,k,a,b,e, alpha, beta, epsilon):
    start = a
    end = b
    #e stands for the error for the approximation
    gap = b-a
    while gap>e:
        
        m = (start+end)/2
        if exact_fun(m,alpha,beta,epsilon)>k:
            end = m
        else:
            start = m
        gap = end-start
    return m

#Now we generate the fine grid

def generate_grid(antecedent, antecedent_k, exact_fun,a,b,e,alpha,beta,epsilon):
    x = np.zeros(len(antecedent))
    for i in range (len(antecedent)):
        x[i] = antecedent_k(exact_fun,antecedent[i],a,b,e,alpha,beta,epsilon)
    return x
  
x = generate_grid(y, antecedent_k, exact_fun,0,1,10e-10,alpha,beta,epsilon)
#print(x)
x[0] = 0
x[-1] = 1


#We construct our matrice A and the vector F

def generate_A_interior(epsilon, m, x):
    A = np.zeros([m+2, m+2])
    for i in range(1,m+1):
        ci = 2/((x[i]-x[i-1])*(x[i]-x[i+1]))
        ci_left = 2/((x[i-1]-x[i])*(x[i-1]-x[i+1]))
        ci_right = 2/((x[i+1]-x[i])*(x[i+1]-x[i-1]))
        A[i, i-1] = epsilon*ci_left + 1/(x[i+1]-x[i-1])
        A[i, i] = epsilon*ci
        A[i, i+1] = epsilon*ci_right - 1/(x[i+1]-x[i-1])
    A[0,0] = 1
    A[-1,-1] = 1
    
    return A
A = generate_A_interior(epsilon, m, x)

def fun(x):
    return -1

def generate_F_interior(m, x, fun, alpha, beta):
    F = np.zeros(m+2)
    for i in range(1,m+1):
        F[i] = fun(x[i])
    F[0] = alpha
    F[-1] = beta
    return F

F = generate_F_interior(m, x, fun, alpha, beta)
def draw_approximation(m, alpha, beta, epsilon, fun):
    y = generate_grid_y(alpha, beta, m)
    x = generate_grid(y, antecedent_k, exact_fun,0,1,10e-10,alpha,beta,epsilon)
    A = generate_A_interior(epsilon, m, x)
    F = generate_F_interior(m, x, fun, alpha, beta)
    U = np.linalg.solve(A, F)
    u = exact_fun(x, alpha, beta, epsilon)
    plt.plot(x,U,'r',linewidth=2)
    plt.plot(x,u,'b',linewidth=2)
    Emax = np.max(np.abs(U-u))
    return Emax
#draw_fun(number, alpha, beta, epsilon)
w = interactive(draw_approximation, alpha=fixed(alpha), beta=fixed(beta), fun=fixed(fun), epsilon=fixed(epsilon),\
                m=widgets.IntSlider(min=1,max=100,step=1, value=1.00, description='$\nnumber:$',))
display(w)
#Let's find the best value of m

Error = np.zeros(100)
for i in range(100):
    Error[i] = draw_approximation(i, alpha, beta, epsilon, fun)

Error = Error[1:len(Error)]
bestError, best_m = min(Error), np.argmin(Error)+1
print(bestError, best_m)<jupyter_output>0.0006843223880734772 99<jupyter_text>Poisson Problem To solve this first question, we use the 5-point Laplacian. We will order unknowns using natural rowwise ordering. The matrix A looks like (Plot that matrix in the final document)<jupyter_code>#Now we compute the matrices T and D

def generate_TD(N, dx, dy):
    T = np.zeros([N,N])
    a = (dy/dx)**2
    b = -2*(1 + a)
    for i in range(N):
        T[i,i] += b
        if (i < N-1):
            T[i,i+1] += a
        if (i > 0):
            T[i,i-1] += a
    D = np.identity(N)
    return T, D

#We can have a loook on the matrix A


def assemble_matrix_A(dx, dy, N, M):
    T, D = generate_TD(N, dx, dy)
    A = np.zeros([N*M, N*M])
    for j in range(M):
        A[j*N:(j+1)*N,j*N:(j+1)*N] += T
        if (j < M-1):
            A[j*N:(j+1)*N,(j+1)*N:(j+2)*N] += D
        if (j > 0):
            A[j*N:(j+1)*N,(j-1)*N:j*N] += D
    return A

N = 3
M = 3
dx = 1./(N+1)
dy = 1./(M+1)
T, D = generate_TD(N, dx, dy)
#print (T)
A = assemble_matrix_A(dx, dy, N, M)
#print (A)
#plt.spy(A)
#plt.show()


#We set boundary conditons
#We assume x and y are in the range [0,1]

def gL(y):
    return np.exp(y/2)

def gR(y):
    return np.exp(1+y/2)

def gB(x):
    return np.exp(x)

def gT(x):
    return np.exp(x+1/2)

#Right side

def f(x,y):
    return 1.25*(np.exp(x+y/2))

#We define the vector of the right side. We must take into consideration the boundary conditions

def assemble_vector_b(x, y, dx, dy, N, M, gL, gR, gB, gT,f):
    b = np.zeros(N*M)
    # Left BCs
    for j in range(M):
        b[(j)*N] += -(dy/dx)**2*gL(y[j+1])
        
    for j in range(M):
        b[(j)*N+N-1] += -(dy/dx)**2*gR(y[j+1]) 
    
    # Bottom BCs
    for i in range(N):
        b[i] += -gB(x[i+1])
    
    # Top BCs:
    for i in range(N):
        b[(M-1)*N+i] += -gT(x[i+1])
        
    for i in range(N):
        for j in range(1,M+1):
            b[(j-1)*N+i] += dy*dy*f(x[i+1],y[j])
            
    return b
    


x = np.linspace(0, 1, N+2)
y = np.linspace(0, 1, M+2)
b = assemble_vector_b(x, y, dx, dy, N, M, gL, gR, gB, gT,f)
v = np.linalg.solve(A,b)
def exact_solution(x,y):
    return np.exp(x+y/2)

#We can compute the vector of the exact solution
u_exc = np.zeros(N*M)
for i in range(N):
    for j in range(1,M+1):
        u_exc[(j-1)*N+i] += exact_solution(x[i+1],y[j])
print(u_exc, v)<jupyter_output>[1.45499141 1.86824596 2.39887529 1.64872127 2.11700002 2.71828183
 1.86824596 2.39887529 3.08021685] [1.4554158  1.86886036 2.39942128 1.64930019 2.11783982 2.71902002
 1.86872606 2.39956934 3.08084177]
{'fraction_non_alphanumeric': 0.161619283065513, 'fraction_numerical': 0.052379480840543884, 'mean_word_length': 3.0180012414649284, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '306358', 'license': 'MIT', 'n_tokens_mistral': 2938, 'n_tokens_neox': 2724, 'n_words': 703, 'path': 'Other/Project1.ipynb', 'repo_name': 'xiaozhouli/Jupyter', 'size': '6472'}