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

import math

import matplotlib.pyplot as plt
import numpy as np



def getK(l):
    k=np.array([[(l**3*(l**2+10))/30,(l**4*(l**2+10))/20],[(l**4*(l**2+10))/20,(4*l**5*(2*l**2+21))/105]])
    return k
def getF(l):
    f=np.array([(l**3/6),(l**4/4)])
    return f
def getN(x,l,n):
    if n==0:
        N=np.array([(x**2-l*x),(x**3-l**2*x)])
    elif n==1:
        N=np.array([(2*x-l),(3*x**2-l**2)])
    return N
def getya(x,l,a,n):
    n_x=len(x)
    ya=np.zeros(n_x)
    for i in range(1,n_x):
        xi=x[i-1]
        if n==0:
            N=getN(xi,l,0)
        elif n==1:
            N=getN(xi,l,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):
        ye[i-1]=math.exp(x[i-1])*a+b*math.exp(-x[i-1])-1
    return ye

l=1
k=getK(l)
f=getF(l)
a=np.matmul(np.linalg.inv(k),f)
x=np.arange(0,1,0.1)
ya=getya(x,l,a,0)
ye=getye(x,l)
plt.plot(x,ya,'r')
plt.plot(x,ye,'b')
plt.show()