Python实现基于变分原理的数值求解：谱方法

import numpy as np
import math

import matplotlib.pyplot as plt


def getK(L,N):
    n=np.arange(1,N+1)
    K=L/2+n**2*math.pi**2/2/L
    K=np.diag(K)
    return K

def getF(L,N):
    n=np.arange(1,N+1).T
    f=np.zeros(N)
    for i in range(0,N+1,2):
        f[i]=(-2)*L/n[i]/math.pi
    return f

def getN(x,L,Nb,n):
    if n==0:
        N=np.sin((math.pi*np.arange(1,Nb+1).T*x)/L)
    elif n==1:
        N=(math.pi*n*np.cos((math.pi*np.arange(1,Nb+1).T*x)/L))/L
    return N
def getya(x,L,a,n,Nb):
    N_x=len(x)
    ya=np.zeros(N_x)
    for i in range(1,N_x+1):
        xi=x[i-1]
        if n==0:
            N=getN(xi,L, Nb,0)
        elif n==1:
            N=getN(xi,L,Nb,1)
        ya[i-1]=np.matmul(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
L=1
N=3
K=getK(L,N)
f=getF(L,N)
a=np.linalg.inv(K)@f
x=np.arange(0,1.1,0.1)
ya=getya(x,L,a,0,N)
ye=getye(x,L)
plt.plot(x,ya,'r')
plt.plot(x,ye,'b')
plt.show()