基函数做最佳平方逼近

使用基函数$1,x,x^2,\dots ,x^n$做最佳平方逼近

原函数$$f(x)=\sqrt{1+x^2}$$
代码如下：

import numpy as np
import math
import matplotlib.pyplot as plt #plt用于显示图片
import sympy
from scipy import integrate
#本程序实现了以1,x,x^2,...x^n为phi的最佳平方逼近
#对应课本例题p86页例5
#scipy.integrate 是积分程序
#scipy.linalg.inv():计算方阵的逆
#sympy是一个符号计算的python库
#对比其他的计算机代数系统，使用sympy需要明确声明符号变量
#求d_i,其中B为对应的phi_k,f为原函数。对f使用phi进行最佳平方逼近

def sum_d(a,b,n):#求d
    x = sympy.symbols('x')
    f = sympy.sqrt(1+pow(x,2))
    B = [sympy.Symbol('m') for w in range(n)]
    v = np.zeros(n)
    for i in range(n):
        B[i]=pow(x,i)*f #修改phi函数
        v[i]=sympy.integrate(B[i],(x,a,b))
    return v



def sum_f(a,b,n):#求H
    y=err=g = np.zeros((n,n,2))
    for i in range(n):
        for j in range(n):
            y[i][j]=integrate.quad(lambda x:pow(x,i+j), a, b)
    #
    return y[...,0]


def can_f(a,b,n):#求a_0,...,a_n
    q,g=np.zeros((n,n))
    q=sum_f(a,b,n)#求H(希尔伯特矩阵)
    g=np.linalg.inv(q)#求逆
    s=v=np.zeros(n)
    s=sum_d(a,b,n)#求d_i
    v=np.dot(g,s)#H^-1*d
    return v

z=can_f(0,1,2)#求a_0,...,a_n
str='%-4.3fx^1+%-4.3f'%(z[1],z[0])#不知道怎样自动生成公式
print(str)

使用勒让德多项式做[-1,1]上的最佳平方逼近

import numpy as np
import math
import matplotlib.pyplot as plt #plt用于显示图片
import sympy 
from scipy import integrate
import Legendr
#本程序实现了以勒让德多项式为phi的最佳平方逼近
#对应课本例题p89页例6
#scipy.integrate 是积分程序
#scipy.linalg.inv():计算方阵的逆
#sympy是一个符号计算的python库
#对比其他的计算机代数系统，使用sympy需要明确声明符号变量
#求d_i,其中B为对应的phi_k,f为原函数。对f使用phi进行最佳平方逼近

x = sympy.symbols('x')#全局变量x为字符常量

def Lf(n):
    B = [sympy.Symbol('m') for w in range(n)]
    for i in range(n):
        B[i]=Legendr.P(i)
    return B


def sum_d(a,b,n):#求d
    f = math.e**x
    B = [sympy.Symbol('m') for w in range(n)]
    v = np.zeros(n)
    for i in range(n):
        B[i]=Lf(n)[i]*f #修改phi函数
        v[i]=sympy.integrate(B[i],(x,a,b))
    return v

#生成矩阵



def can_f(a,b,n):#求a_0,...,a_n
    q=np.zeros((n,n))
    for i in range(n):
        for j in range(n):
            if (i==j):
                q[i][j]=(2*i+1)/2
            else:
                q[i][j]=0
    print(q)
    s=v=np.zeros(n)
    s=sum_d(a,b,n)#求d_i
    v=np.dot(q,s)#H^-1*d
    return v
#求三次最佳平方逼近多项式，n应该取4

z=can_f(-1,1,4)
t=sympy.simplify(z[0]*Lf(4)[0]+z[1]*Lf(4)[1]+z[2]*Lf(4)[2]+z[3]*Lf(4)[3])#代入多项式并简化
print(t)

#赋值计算
print(math.exp(1))
print(t.subs(x,1))

其中，Legendr.py为

import sympy

x = sympy.symbols('x')

def P(point):
    if point == 0:
        return 1
    elif point == 1:
        return x
    p0 = 1
    p1 = x
    for i in range(point - 1):
        temp = ((2 * i + 3) * x * p1 - (i + 1) * p0) / (i + 2)
        p0 = p1
        p1 = temp
    return p1

使用勒让德多项式做[a,b]上的最佳平方逼近并绘图

代码如下：

import numpy as np
import math
import matplotlib.pyplot as plt #plt用于显示图片
import sympy 
from scipy import integrate
import Legendr
import matplotlib.pyplot as plt
#本程序实现了以勒让德多项式为phi的最佳平方逼近
#对应课本例题p89页例6
#scipy.integrate 是积分程序
#scipy.linalg.inv():计算方阵的逆
#sympy是一个符号计算的python库
#对比其他的计算机代数系统，使用sympy需要明确声明符号变量
#求d_i,其中B为对应的phi_k,f为原函数。对f使用phi进行最佳平方逼近

x = sympy.symbols('x')#全局变量x为字符常量
g = sympy.symbols('g')

def Lf(n):
    B = [sympy.Symbol('m') for w in range(n)]
    for i in range(n):
        B[i]=Legendr.P(i)
    return B


def sum_d(a,b,n):#求d
    #与上一程序相比变化之处(做了变量替换)
    f = math.e**((b-a)*x/2+(b+a)/2)
    B = [sympy.Symbol('m') for w in range(n)]
    v = np.zeros(n)
    for i in range(n):
        B[i]=Lf(n)[i]*f #修改phi函数
        v[i]=sympy.integrate(B[i],(x,-1,1))
    return v

#生成矩阵

def can_f(a,b,n):#求a_0,...,a_n
    q=np.zeros((n,n))
    for i in range(n):
        for j in range(n):
            if (i==j):
                q[i][j]=(2*i+1)/2
            else:
                q[i][j]=0
    s=v=np.zeros(n)
    s=sum_d(a,b,n)#求d_i
    #print(s)
    v=np.dot(q,s)#H^-1*d
    return v

def ttox(r,a,b):
    s=r.subs(x,(2*x-a-b)/(b-a))
    return s

#求三次最佳平方逼近多项式，n应该取4
def valuey(a,b,n,c):
    z=can_f(a,b,n)
    t=sympy.simplify(z[0]*Lf(4)[0]+z[1]*Lf(4)[1]+z[2]*Lf(4)[2]+z[3]*Lf(4)[3])#代入多项式并简化
    s=ttox(t,1,5).subs(x,c)
    return s

#画图像
def drawv(a,b,n):
    s=np.arange(a,b,0.01)
    vy=[math.exp(s[i]) for i in range(len(s))]
    y=[valuey(a,b,n,s[i]) for i in range(len(s))]
    plt.plot(s,y,linewidth=3, color='lawngreen')
    plt.plot(s,vy,linewidth=3, color='cyan')
    plt.show()

drawv(1,5,4)