Python实现基于变分原理的数值求解：无网格法

import math

import matplotlib.pyplot as plt
import numpy as np
L=1
dx=0.1
X=np.arange(0,L+dx,dx)
h=1.1

X_M=np.arange(0,1+dx,dx)

def getK(X,h,X_M):
    N_M=len(X_M)
    Nd=len(X)
    K=0
    L=max(X)-min(X)
    
    for i in range(1,N_M+1):
        xi=X_M[i-1]
        Nx=getN(xi,X,h,1)
        N=getN(xi,X,h,0)
        K=K+(Nx@Nx.T+N@N.T)
    
    K=K*(L/N_M)
    return K

def getf(X,h,X_M):
    N_M=len(X_M)
    Nd=len(X)
    f=0
    L=max(X)-min(X)
    
    for i in range(1,N_M+1):
        xi=X_M[i-1]
        N=getN(xi,X,h,0)
        f=f+N
        
    f=-f*(L/N_M)
    return f

def getNn(x,h,n):
    if abs(x)>3*h:
        Nn=0
    else:
        if n==0:
            Nn=math.exp(-x**2/h**2)
        elif n==1:
            Nn=-(2*x*math.exp(-x**2/h**2))/h**2
    
    return Nn

def getN(x,X,h,n):
    Nd=len(X)
    N=np.zeros(Nd)
    for i in range(1,Nd+1):
        xn=X[i-1]
        N[i-1]=getNn(x-xn,h,n)
    return N

def getya(x,X,h,a,n):
    N_x=len(x)
    ya=np.zeros(N_x)
    for i in range(1,N_x+1):
        xi=x[i-1]
        N=getN(xi,X,h,n)
        ya[i-1]=N.T@a
    return ya

def getye(x,l):
    a=(1-math.exp(-l))/(math.exp(l)-math.exp(-l))
    b=(math.exp(l)-1)/(math.exp(l)-math.exp(-l))
    n_x = len(x)
    ye = np.zeros(n_x)
    for i in range(1, n_x+1):
        ye[i-1]=math.exp(x[i-1])*a+b*math.exp(-x[i-1])-1
    return ye

K=getK(X,h,X_M)
f=getf(X,h,X_M)
alfa=10000
N0=getN(0,X,h,0)
NL=getN(L,X,h,0)
K=K+alfa*N0@N0.T+alfa*NL@NL.T
a=K/f
x=np.arange(0,1.1,0.1)
ya=getya(x,X,h,a,0)
ye=getye(x,L)
# plt.plot(X,ya,'r-o')
plt.plot(X,ye,'b')
plt.show()