vandermonda方程的求解

一、实验目的
求解线性方程组 ，其中V是vandermonda行列式。
二、实验过程
1.遇到的问题：
如果按照一般线性方程组的求解方法或一般矩阵的求解逆方法进行vandermonda矩阵的有关计算，所需运算量为O(n^3)。这样大的运算量显示这种方法没有充分利用该类矩阵的结构特点与性质。另外，vandermonda矩阵的条件数一般都很大，从而相应的求解问题是病态的。
2.采用的解决方法：
研究在某些条件下稳定性好的算法，以求得高精度的解。采用Bjorck - Pereyra 算法求解通常的vandermonda 方程。
3.Bjorck - Pereyra 算法步骤：（参考文献：《范德蒙矩阵类的快速算法》徐仲，西北工业大学出版社，第二章的第2小节）

#Bjorck - Pereyra 算法（vandermonda方程的快速算法）
#即求解 Vx=b
import numpy as np

#step1:生成vandermonda方程
a = input('请输入ai用","分隔:')
a = list(map(int, a.split(",")))     
#input接受的是字符串，先用.split将其进行分割，再用map函数将int函数迭代作用到每个元素上，最后用list使其成为一个整型数组
a = np.array(a)
n = len(a)
V = np.vander(a)         #V是方程组的系数矩阵，其行列式为vandermonda行列式
V = np.transpose(V)
V = V[::-1]              #使用.vander()函数得到的数组不是常规的、我想要的数组，必须要经过调整（unbeliveable）
b = input('请输入bi用","分隔且须与ai一样长:')
b = list(map(int, b.split(",")))
b[0:n]     #截取前n个值，确保最后得到的是一个方程组
reb = np.array(b)
print("所输入的vandermonda矩阵（系数矩阵V）为：\n",V)
print("所输入的方程组的解b为：",b)

#step2:求dij
d = np.zeros((n,n))
for j in range(0,n):                #对应算法step2,line1
    d[j][0] = b[j]     
for k in range(0,n-1):        #对应算法step2,line2
    for j in range(k, n):
        d[j][k+1] = d[j][k] - a[k]*d[j-1][k]

#step3:求解y
z = np.zeros((n,n))   #用于存放中间变量
y = np.zeros((n,n))
for j in range(0,n):
    y[j][n-1] = d[j][j]
#add by myself (in order to the problem of yij)
for j in range(n):
    for i in range(j):
        y[i][j] = y[i][n-1]   #这两个循环在于补全yij,否则在下面推yj1的部分会由于数据缺失导致部分结果不准确
    
for k in range(n-2,-1,-1):
    z[k][k+1] = y[k][k+1]
    for j in range(k+1, n):
        z[j][k+1] = y[j][k+1]/(a[j] - a[j-k-1])  #注意除的第二部分由a[j-k]改成了a[j-k-1],
                                           #原本的算法k是从1开始的，而编程习惯是以0开始，从而会导致分母为0的情况。
    for j in range(k, n-1):
        y[j][k] = z[j][k+1] - z[j+1][k+1]
    y[n-1][k] = z[n-1][k+1]

#输出结果
print("所求方程组的解y为：",y[:,0])
print("验算部分V*y:",np.dot(V,y[:,0])) 

结果：