插值法

实现范德蒙德多项式插值、拉格朗日插值、牛顿插值、分段线性、分段三次 Hermite 插值，并完成各方法之间的对比。
输入：插值区间 [a,b][a, b]，参数 c,d,e,fc, d, e, f 作为标准函数 f(x)=c⋅sin⁡dx+e⋅cos⁡fxf(x) = c \cdot \sin dx + e \cdot \cos fx 的值，参数 n+1n+1 作为采样点的个数，参数 mm 作为实验点的个数。
要求：在区间 [a,b][a, b] 上均匀采集个采集点，利用这 n+1n+1 个采集点，分别使用范德蒙德多项式插值、拉格朗日插值、牛顿插值、分段线性、分段三次 Hermite 插值进行插值，求出 L(x)L(x)，之后再选取 mm 个点作为实验点，计算在这 mm 个实验点上插值函数 L(x)L(x) 与目标函数 f(x)f(x) 的平均误差。同时对比各插值方法之间的精度差异。
输出：对比函数曲线，平均误差。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
import numpy as np
import matplotlib.pyplot as plt

# 范德蒙德多项式插值
def vandermonde_interpolation(x, y, x_values):
    n = len(x)
    coefficients = []
    c = []
    for i in range(n):
        coeff = []
        for j in range(n):
            coeff.append(x[i]**j)
        coefficients.append(coeff)
        c.append(y[i])
    A = np.array(coefficients)
    inv_A = np.linalg.inv(A)
    a = inv_A.dot(c)
    sum = 0
    for i in range(n):
        sum += a[i] * x_values ** i
    return sum

# 拉格朗日插值
def lagrange_interpolation(x, y, x_values):
    n = len(x)
    result = 0
    for i in range(n):
        p = y[i]
        for j in range(n):
            if j != i:
                p *= (x_values - x[j]) / (x[i] - x[j])
        result += p
    return result

# 牛顿插值
def newton_interpolation(x, y, x_values):
    n = len(x)
    coefficients = []
    for i in range(n-1, 0, -1):
        divided_diff = (y[i] - y[i-1]) / (x[i] - x[i-1])
        coefficients.append(divided_diff)
        for j in range(n-i-1, 0, -1):
            divided_diff = (-divided_diff + coefficients[j-1]) / (x[n-j] - x[i-1])
            coefficients[j-1] = divided_diff
    coefficients.append(y[0])
    result = coefficients[0]
    for i in range(n-2, -1, -1):
        result = result * (x_values - x[i]) + coefficients[n-1-i]
    return result

# 分段线性插值
def piecewise_linear_interpolation(x, y, x_values):
    n = len(x)
    result = 0
    for i in range(n-1):
        mask = (x[i] <= x_values) & (x_values < x[i+1])
        slope = (y[i+1] - y[i]) / (x[i+1] - x[i])
        result += ((y[i] + slope * (x_values - x[i]))*mask)
    if (x_values == x[n-1]):
        result = y[n-1]
    return result

# 分段三次Hermite插值
def piecewise_cubic_hermite_interpolation(x, y, yy, x_values):
    n = len(x)
    result = 0
    for i in range(n-1):
        if ((x[i] <= x_values) & (x_values < x[i+1])):
            ai = (1+2*(x_values-x[i])/(x[i+1]-x[i]))*((x_values-x[i+1])/(x[i]-x[i+1]))**2
            bi = (x_values - x[i])*((x_values-x[i+1])/(x[i]-x[i+1]))**2
            ai1 = (1+2*(x_values-x[i+1])/(x[i]-x[i+1]))*((x_values-x[i])/(x[i+1]-x[i]))**2
            bi1 = (x_values - x[i+1])*((x_values-x[i])/(x[i+1]-x[i]))**2
            result += (y[i] * ai +yy[i] * bi +y[i+1] * ai1 +yy[i+1] * bi1)
    if (x_values == x[n-1]):
        result = y[n-1]
    return result

# 目标函数
def target_function(c,d,e,f,x):
    return [(c*np.sin(d*val)+e*np.cos(f*val)) for val in x] 

# 目标函数导函数
def derivative_function(c,d,e,f,x):
    return [(c*d*np.cos(d*val)-e*f*np.sin(f*val)) for val in x]

# 计算平均误差
def compute_average_error(f, g, x_values):
    return sum([abs(f[i] - g[i]) for i in range(len(x_values))]) / len(x_values)

# 设置插值区间和参数
a = float(input())
b = float(input())
c = float(input())
d = float(input())
e = float(input())
f = float(input())
num_samples = int(input())
num_experiments = int(input())
x_values = [a + (b - a) * i / (num_samples - 1) for i in range(num_samples)]
ys = target_function(c,d,e,f,x_values)
yy = derivative_function(c,d,e,f,x_values)

# 计算实验点
experiment_points = [a + (b - a) * i / (num_experiments - 1) for i in range(num_experiments)]

# 计算各插值方法的插值结果和平均误差
vandermonde_interpolated = [vandermonde_interpolation(x_values, ys, val) for val in experiment_points]
vandermonde_error = compute_average_error(target_function(c,d,e,f,experiment_points), vandermonde_interpolated, experiment_points)

lagrange_interpolated = [lagrange_interpolation(x_values, ys, val) for val in experiment_points]
lagrange_error = compute_average_error(target_function(c,d,e,f,experiment_points), lagrange_interpolated, experiment_points)

newton_interpolated = [newton_interpolation(x_values, ys, val) for val in experiment_points]
newton_error = compute_average_error(target_function(c,d,e,f,experiment_points), newton_interpolated, experiment_points)

piecewise_linear_interpolated = [piecewise_linear_interpolation(x_values, ys, val) for val in experiment_points]
piecewise_linear_error = compute_average_error(target_function(c,d,e,f,experiment_points), piecewise_linear_interpolated, experiment_points)

piecewise_cubic_hermite_interpolated = [piecewise_cubic_hermite_interpolation(x_values, ys, yy, val) for val in experiment_points]
piecewise_cubic_hermite_error = compute_average_error(target_function(c,d,e,f,experiment_points), piecewise_cubic_hermite_interpolated, experiment_points)

fig = plt.figure(num = 1,dpi = 120)
ax = plt.subplot(1,1,1)
 # 坐标轴
ax = plt.gca()  # get current axis 获得坐标轴对象
ax.spines['right'].set_color('none')  # 将右边 边沿线颜色设置为空 其实就相当于抹掉这条边
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
# 设置中心的为（0，0）的坐标轴
ax.spines['bottom'].set_position(('data', 0))  # 指定 data 设置的bottom(也就是指定的x轴)绑定到y轴的0这个点上
ax.spines['left'].set_position(('data', 0))

x=list(np.arange(a,b,0.01))#此处可调整自变量取值范围，以便选择合适的观察尺度
y=[]
y1=[]
y2=[]
y3=[]
y4=[]
y5=[]
for i in range(len(x)):
    y = target_function(c,d,e,f,x)
    y1.append(vandermonde_interpolation(x_values, ys, x[i]))
    y2.append(lagrange_interpolation(x_values, ys, x[i]))
    y3.append(newton_interpolation(x_values, ys, x[i]))
    y4.append(piecewise_linear_interpolation(x_values, ys, x[i]))
    y5.append(piecewise_cubic_hermite_interpolation(x_values, ys, yy, x[i]))
    

ax.plot(x,y,label = "Target Function",color ="blueviolet")
ax.plot(x_values,ys, marker = "*",linestyle = "", color = "blueviolet")
ax.plot(x,y1,label = "vandermonde interpolation\n average error=%f"%vandermonde_error,color ="red")
ax.plot(experiment_points, vandermonde_interpolated, marker = "o",linestyle = "", color = "red")
ax.plot(x,y2,label = "lagrange interpolation\n average error=%f"%lagrange_error,color ="yellow")
ax.plot(experiment_points, lagrange_interpolated, marker = "o",linestyle = "", color = "yellow")
ax.plot(x,y3,label = "newton interpolation\n average error=%f"%newton_error,color ="green")
ax.plot(experiment_points, newton_interpolated, marker = "o",linestyle = "", color = "green")
ax.plot(x,y4,label = "piecewise linear interpolation\n average error=%f"%piecewise_linear_error,color ="blue")
ax.plot(experiment_points, piecewise_linear_interpolated, marker = "o",linestyle = "", color = "blue")
ax.plot(x,y5,label = "piecewise cubic hermite interpolation\n average error=%f"%piecewise_cubic_hermite_error,color ="purple")
ax.plot(experiment_points, piecewise_cubic_hermite_interpolated, marker = "o",linestyle = "", color = "purple")
#ax.set_xlim(0,2)
#plt.draw()
plt.legend()
plt.show()

函数逼近

实现最佳平方逼近与最小二乘拟合，并完成两种方法之间的对比。
输入：函数区间 [a,b][a, b]，参数 cc 作为标准函数f(x)=11+cx2f(x)=\frac{1}{1+cx^2}的值，参数 kk 作为所构造的逼近多项式的次数 (k=1,2,3)(k=1,2,3) 。参数 n+1n+1 作为采样点的个数，参数 mm 作为实验点的个数。
要求：要求选用勒让德正交多项式作最佳平方逼近；在区间 [a,b][a, b] 上均匀采集 nn 个采集点，利用这 n+1n+1 个采集点，计算采集点上的函数值，构造最小二乘拟合多项式函数。之后再选取 mm 个点作为实验点，计算在这 mm 个实验点上所构造的逼近函数与给定的目标函数 f(x)f(x) 的平均误差。同时对比两种逼近方法之间的精度差异。
输出：对比函数曲线，平均误差。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
from sympy import *
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
def qpow(a,n):
    ans=1
    while n!=0:
        if n&1:
            ans*=a
        a*=a
        n>>=1
    return ans
def F(x):
    res = 1/(1+c*x*x)
    return res
def cntpf():
    t = symbols('t')
    # (b-a)*t/2+(b+a)/2
    fp[0] = integrate(1 / (1 + c * ((b - a) * t / 2 + (b + a) / 2) * ((b - a) * t / 2 + (b + a) / 2)),(t, -1, 1)).evalf()
    fp[1] = integrate((1 / (1 + c * ((b - a) * t / 2 + (b + a) / 2) * ((b - a) * t / 2 + (b + a) / 2))) * (t),(t, -1, 1)).evalf()
    fp[2] = integrate((1 / (1 + c * ((b - a) * t / 2 + (b + a) / 2) * ((b - a) * t / 2 + (b + a) / 2))) * (3 * t * t / 2 - 1 / 2),(t, -1, 1)).evalf()
    fp[3] = integrate((1 / (1 + c * ((b - a) * t / 2 + (b + a) / 2) * ((b - a) * t / 2 + (b + a) / 2))) * (5 * t * t * t / 2 - 3 * t / 2),(t, -1, 1)).evalf()
def bsa(inp):
    res=0
    inp =(2 * inp -a - b) / (b - a)
    for i in range(0,k):
        res+=(2*i+1)*fp[i]*qpow(inp,i)/2
    return res
def setup():
    for i in range(0,100):
        vispa[i]=0
        vispb[i]=0
def initx():
    d=(b-a)/n
    x[0]=a
    for i in range(1,n):
        x[i]=x[i-1]+d
def compa(i):
    if vispa[i]==1:
        return  pa[i]
    up=0
    down=0
    for j in range(0,n):
        up+=comP(i-1,x[j])*comP(i-1,x[j])*x[j]
        down += comP(i - 1, x[j]) * comP(i - 1, x[j])
    res=up/down
    vispa[i] = 1
    pa[i] = res
    return res
def compb(i):
    if vispb[i]==1:
        return pb[i]
    up = 0
    down = 0
    for j in range(0, n):
        up += comP(i, x[j]) * comP(i, x[j])
        down += comP(i - 1, x[j]) * comP(i - 1, x[j])
    res=up/down
    vispb[i] = 1
    pb[i] = res
    return res
def comP(i,inp):
    res=0
    if i==0:
        res=1
    elif i==1:
        res = (inp -compa(1))*comP(0, inp)
    else:
        res = (inp -compa(i))*comP(i - 1, inp)-compb(i - 1) * comP(i - 2, inp)

    return res
def lsf(inp):
    initx()
    res=0
    for i in range(0,k):
        up=0
        down=0
        for j in range(0,n):
            up += F(x[j]) * comP(i, x[j])
            down += comP(i, x[j]) * comP(i, x[j])
        res += comP(i, inp)*up / down
    return res
def bsaerror(inp):
    return abs(bsa(inp)-F(inp))
def lsferror(inp):
    return abs(lsf(inp)-F(inp))

#全局变量
x=[None]*200
fp=[None]*50
vispa=[0]*100
vispb=[0]*100
pa=[None]*200
pb=[None]*200

#主函数段
a, b, c = map(int, input("输入a,b,c: ").split())
n , k = map(int, input("输入n,k: ").split())
m=input("请输入测试组数：")
tst = input("请输入测试用例：").split()

setup() #初始化fp数组
cntpf() #初始化fp数组

bsax=range(int(a),int(b))
bsay=[]
bsae=0
for i in bsax:
    bsay.append(bsa(i))
    
lsfx=range(int(a),int(b))
lsfy=[]
lsfe=0
for i in lsfx:
    lsfy.append(lsf(i))
    
nx=range(int(a),int(b))
ny=[]
for i in nx:
    ny.append(F(i))
    
for i in tst:
    bsae+=bsaerror(float(i))
    lsfe+=lsferror(float(i))
bsae/=float(m)
lsfe/=float(m)
print('最佳平方逼近的平均误差为%f'%bsae)
print('最小二乘拟合的平均误差为%f'%lsfe)
plt.plot(bsax, bsay, label='最佳平方逼近', color='red')
plt.plot(lsfx, lsfy, label='最小二乘拟合', color='blue')
plt.plot(nx, ny, label='原函数', color='orange')
plt.legend()
plt.title('函数逼近对比函数')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

数值积分

实现复化梯形公式和龙贝格算法计算积分，并完成两种方法之间的精度对比。
输入：函数区间 [a,b][a, b]，被积函数为 f(x)=xlnxf(x) =\sqrt{x}lnx，参数 hh 作为步长。参数 εε 作为要求满足的精度条件。
要求：取不同的步长 hh，要求用复化梯形公式和龙贝格算法分别计算积分值计算当精度达到 εε 时，所需要等分积分区间的次数（假设每次都是二等分）及 hh 的大小。当达到精度要求时，对比两种方法需要划分次数及步长 hh 的大小。
输出：数值积分计算结果，划分次数，步长 hh 的大小。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
import math
#初始化变量
vis = [[0 for i in range(1000)] for j in range(1000)]
T = [[0 for i in range(1000)] for j in range(1000)]

def f(x):
    res=math.sqrt(x)*math.log(x,math.e)
    return res

def F(x):
    x3=x*x*x
    res=2*math.log(x,math.e)*math.sqrt(x3)/3-4*math.sqrt(x3)/9
    return res

#使用牛顿莱布尼茨公式直接积分
def Fitg(a,b):
    res=F(b)-F(a)
    return  res

#使用复合梯形求积公式直接积分
def ctqf(n):
    h=(b-a)/n
    res=0
    for k in range(n):
        res+=f(a+k*h)+f(a+(k+1)*h)
    res=res*h/2
    return res

#计算T的值
def cntT(m,k):
    if(vis[m][k]==1):
        return T[m][k]
    else:
        vis[m][k] == 1
        if(m==0 and k==0):
            T[m][k]=(b-a)*(f(a)+f(b))/2
        elif(m==0):
            n=1<<(k-1)
            h=(b-a)/n
            sm=0
            for i in range(n):
                sm+=f(a+(2*i+1)*h/2)
            sm=sm*h/2
            T[m][k]=cntT(m,k-1)/2+sm
        else:
            fourm=math.pow(4,m)
            T[m][k]=fourm*cntT(m-1,k+1)/(fourm-1)-cntT(m-1,k)/(fourm-1)
    return T[m][k]

#使用龙贝格积分求积
def Romberg(m):
    return cntT(m,0)

a,b =map(float,input("请输入积分区间a,b: ").split())
e=float(input("请输入精度条件e: "))

for k in range(1000):
    n=1<<k
    eps=abs(ctqf(n)-Fitg(a,b))
    if(eps<e):
        h=(b-a)/n
        print('复合梯形求积公式的计算结果为：%f'%(ctqf(n)))
        print('划分次数为：%d'%(k))
        print('步长为：%f' % (h))
        break

for k in range(1,1000):
    n=1<<k
    eps=abs(Romberg(k)-Romberg(k-1))
    if(eps<e):
        h=(b-a)/n
        print('龙贝格求积公式的计算结果为：%f'%(Romberg(k)))
        print('划分次数为：%d'%(k))
        print('步长为：%f' % (h))
        break

线性方程组的直接解法

1. 编写列主元消元法的通用程序。
输入：矩阵 AA 和向量 bb。
输出：消元后的增广矩阵及方程 Ax=bAx=b 的根值 x∗x^*。
要求1：选取进行测试，打印出上述程序的输出。
要求2：随机生成 nn 阶（nn>=20，具体值自定）方阵 AA，n∗1n*1 维非零向量 bb，求解 xx。

2. 编写使用LU分解法求解线性方程组的通用程序。
输入：矩阵 AA 和向量 bb。
输出：对矩阵 AA 进行 LULU 分解后的 LL 和 UU，以及方程组的根值 x∗x^*。
要求1：选取进行测试，打印输出结果。
要求2：随机生成 nn 阶（nn>=20，具体值自定）方阵 AA，n∗1n*1 维非零向量 bb，求解 xx。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
#include<iostream>
#include<algorithm>
#include<ctime>
#include <random>

using namespace std;
int n;
const int maxn=200;
double inA[maxn][maxn];
double A[maxn][maxn];
double L[maxn][maxn];
double U[maxn][maxn];
double ans[maxn];
double ansy[maxn];
void Pivot();//列主元消元法
void copyA();//将输入矩阵复制给A
void printA();//打印矩阵A
void printb();//打印向量b
void exc(int i,int j);//交换行
void Minus(int i,int j,double mult);//第j行减第i行
void printans();//打印答案
void LU();//LU分解法
void iscorrect();//回代判断答案是否正确
double gerandom(double a,double b);//生成a到bdouble类型的随机数
int main(){
    int T;
     cout<<"please input T :";
     cin>>T;
     while (T--)
     {
         cout << "please input n :";
         cin >> n;
         cout << "enter 1 if you want to use your own input or 2 if you want your input to be random"<<endl;
         int cho;
         cin>>cho;
         if (cho == 1)
         {
             cout << "please input matrix A :" << endl;
             for (int i = 1; i <= n; i++)
             {
                 for (int j = 1; j <= n; j++)
                 {
                     cin >> inA[i][j];
                 }
             }
             cout << "please input vector b :" << endl;
             for (int i = 1; i <= n; i++)
             {
                 cin >> inA[i][n + 1];
             }
         }
         else{
            double a,b;
            cout << "please input the range of the random number :" << endl;
            cin>>a>>b;
            for (int i = 1; i <= n; i++)
             {
                 for (int j = 1; j <= n; j++)
                 {
                    inA[i][j]=gerandom(a,b);
                 }
             }
             for (int i = 1; i <= n; i++)
             {
                inA[i][n + 1]=gerandom(a,b);
             }
         }
         cout << "enter 1-Pivot or 2-LU:" << endl;
         cin>>cho;
         if(cho==1){
            Pivot();
         }
         else{
            LU();
         }
         printans();
         iscorrect();
     }
    return 0;
}
void printA(){
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            cout<<inA[i][j]<<" ";
        }
        cout<<endl;
    }
}
void copyA(){
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n+1;j++){
            A[i][j]=inA[i][j];
        }
    }
}
void exc(int i,int j){
    for(int k=1;k<=n+1;k++){
        swap(A[i][k],A[j][k]);
    }
}
void Minus(int i,int j,double mult){
    for(int k=1;k<=n+1;k++){
        A[j][k]-=A[i][k]*mult;
    }
}
void printans(){
    for(int i=1;i<=n;i++){
        cout<<ans[i]<<" ";
    }
    cout<<endl;
}
void Pivot(){
    copyA();
    for(int t=1;t<=n-1;t++){
        double nman=0;
        int maxi=0;
        for(int i=t;i<=n;i++){
            double tp=abs(A[i][t]);
            if(tp>=nman){
                maxi=i;
                nman=tp;
            }
        }
        exc(t,maxi);
        for(int i=t+1;i<=n;i++){
            double multi=A[i][t]/A[t][t];
            Minus(t,i,multi);
        }
    }
    for(int t=n;t>=1;t--){
        double right=A[t][n+1];
        for(int j=n;j>t;j--){
            right-=A[t][j]*ans[j];
        }
        ans[t]=right/A[t][t];
    }
}
void iscorrect(){
    int tag=1;
    for(int i=1;i<=n;i++){
        double myans=0;
        for(int j=1;j<=n;j++){
            myans+=inA[i][j]*ans[j];
        }
        if(abs(myans-inA[i][n+1])>=1e-4){
            tag=0;
        }
    }
    if(tag){
        cout<<"answer is correct!"<<endl;
    }
    else{
        cout<<"answer is wrong!"<<endl;
    }
}
void LU(){
    copyA();
    for(int i=1;i<=n;i++){
        U[1][i]=A[1][i];
        L[i][i]=1;
    }
    for(int i=2;i<=n;i++){
        L[i][1]=A[i][1]/U[1][1];
    }
    for(int r=2;r<=n;r++){
        for(int i=r;i<=n;i++){
            double tplu=0;
            for(int k=1;k<=r-1;k++){
                tplu+=L[r][k]*U[k][i];
            }
            U[r][i]=A[r][i]-tplu;
        }
        if(r!=n){
            for(int i=r+1;i<=n;i++){
                double tplu=0;
                for(int k=1;k<=r-1;k++){
                    tplu+=L[i][k]*U[k][r];
                }
                L[i][r]=(A[i][r]-tplu)/U[r][r];
            }
        }
    }
    ansy[1]=A[1][n+1];
    for(int i=2;i<=n;i++){
        double tply=0;
        for(int k=1;k<=i-1;k++){
            tply+=L[i][k]*ansy[k];
        }
        ansy[i]=A[i][n+1]-tply;
    }
    ans[n]=ansy[n]/U[n][n];
    for(int i=n-1;i>=1;i--){
        double tpux=0;
        for(int k=i+1;k<=n;k++){
            tpux+=U[i][k]*ans[k];
        }
        ans[i]=(ansy[i]-tpux)/U[i][i];
    }
}
double gerandom(double a,double b){
    static random_device rd;
    static mt19937 gen(rd());
    uniform_real_distribution<double> dis(a, b);
    return dis(gen);
}

线性方程组的迭代解法

1. 编写高斯-塞德尔迭代和SOR迭代的通用程序。
输入：矩阵A和向量b，迭代初值x0，迭代最大步数K，误差控制。对于超松弛迭代，还需输入松弛因子。
输出：迭代步数及方程 Ax=bAx=b 的根值 x∗x^*。
要求：
(1)选取进行测试，取初值x(0)=0，误差控制=10-8，打印出两种迭代方法的输出。
(2)取松弛因子 =i50=\frac{i}{50}，ii=1,2,…99,打印迭代步数，并给出一个最佳值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
const int maxn = 1000;
double a[maxn][maxn];
double b[maxn];
double xk[maxn];
double xk1[maxn];
int n;
int K;
double w;
double eps;
double mf;
double nstep;
void setxk0();//初始化xk0
void movk();//将xk1赋值给xk
double cntloss();//计算detx
void SOR(double w);//超松弛迭代法
void printx();//打印x
void iscorrect();//回代验证答案是否正确
int main(){
    cout<<"input max step K : ";
    cin>>K;
    cout<<"input eps  : ";
    cin>>eps;
    cout<<"input the size of A  : ";
    cin>>n;
    cout<<"matrix A :"<<endl;
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            cin>>a[i][j];
        }
    }
    cout<<"vector b :"<<endl;
    for(int i=1;i<=n;i++){
        cin>>b[i];
    }
    cout<<"enter 1 if you want to test one w"<<endl;
    cout<<"or 2 if you want to find the best w"<<endl;
    cout<<"enter : ";
    cin>>mf;
    if(mf==1){
        cout<<"input Relaxation factor w  : ";
        cin>>w;
        SOR(w);
        printx();
        iscorrect();
    }
    else{
        int bst=K+1;
        double bsw=0;
        for(int i=1;i<=99;i++){
            double ii = (double)i;
            double nw=ii/(double)50;
            cout<<"now i is "<<i<<endl;
            SOR(nw);
            if(nstep<bst){
                bst = nstep;
                bsw=nw;
            }
        }
        cout<<"the best w is "<<bsw<<endl;
    }
    
    return 0;
}
void setxk0(){
    for(int i=1;i<=n;i++){
        xk[i]=0;
    }
}
void movk(){
    for(int i=1;i<=n;i++){
        xk[i]=xk1[i];
    }
}
double cntloss(){
    double res = 0;
    for(int i=1;i<=n;i++){
        double tp = abs(xk1[i]-xk[i]);
        if(tp>res){
            res = tp;
        }
    }
    return res;
}
void printx(){
    for(int i=1;i<=n;i++){
        cout<<xk1[i]<<" ";
    }
    cout<<endl;
}
void iscorrect(){
    int flag = 1;
    for(int i=1;i<=n;i++){
        double s=0;
        for(int j=1;j<=n;j++){
            s+=a[i][j]*xk1[j];
        }
        if(abs(s-b[i])>1e-3){
            flag=0;
        }
    }
    if(flag){
        cout<<"answer is correct!"<<endl;
    }
    else{
        cout<<"answer is wrong!"<<endl;
    }
}
void SOR(double w){
    setxk0();
    for(int step=1;step<=K;step++){
        for(int i=1;i<=n;i++){
            double s1 = 0;
            double s2 = 0;
            for(int j=1;j<=i-1;j++){
                s1+=a[i][j]*xk1[j];
            }
            for(int j=i;j<=n;j++){
                s2+=a[i][j]*xk[j];
            }
            xk1[i]=xk[i]+w*(b[i]-s1-s2)/a[i][i];
        }
        if(cntloss()<eps){
            if(mf==2){
                nstep = step;
                cout<<"the step is "<<step<<endl;
            }
            break;
        }
        else{
            movk();
        }
    }
}

非线性方程（组）的数值解法

1. 编写不动点迭代、斯特芬森加速迭代和牛顿迭代的通用程序。
要求：
(1)设计一种不动点迭代格式，求解函数 f(x)=x2−3x+2−exf(x)=x^2-3x+2-e^x 和 g(x)=x3+2x2+10x−20g(x)=x^3+2x^2+10x-20 的根，要求该迭代格式收敛。然后再使用斯特芬森加速迭代，计算到 |xk-xk-1|<10-8 为止。
(2)用牛顿迭代，同样计算到 |xk-xk-1|<10-8 。输出迭代初值、迭代次数及各次迭代值，比较方法优劣。

1

2. 本章计算实习题3。

1