2023年数值计算方法大作业.doc

目录
第一章 非线性方程求根····················3
1.1迭代法······································3
1.2牛顿法······································4
1.3弦截法······································5
1.4二分法······································6
第二章 插值·······························7
2.1线性插值····································7
2.2二次插值····································8
2.3拉格朗日插值································9
2.4分段线性插值·······························10
2.5分段二次插值·······························11
第三章 数值积分··························13
3.1复化矩形积分法·····························13
3.2复化梯形积分法·····························14
3.3辛普森积分法·······························15
3.4变步长梯形积分法···························16
第四章 线性方程组数值法··················17
4.1约当消去法·································17
4.2高斯消去法·································18
4.3三角分解法·································20

4.4雅可比迭代法································21
4.5高斯—赛德尔迭代法··························23
第五章 常积分方程数值法··················25
5.1显示欧拉公式法······························25
5.2欧拉公式预测校正法··························26
5.3改善欧拉公式法······························27
5.4四阶龙格—库塔法····························28

数值计算措施
第一章 非线性方程求根

1.1迭代法
程序代码：
Private Sub Command1_Click()
x0 = Val(InputBox("请输入初始值x0"))
ep = Val(InputBox(请输入误差限ep))
f = 0
While f = 0
X1 = (Exp(2 * x0) - x0) / 5
If Abs(X1 - x0) < ep Then
Print X1
f = 1
Else
x0 = X1
End If
Wend
End Sub

例：求f(x)=e2x-6x=0在x=0.5附近旳根（ep=10-10）

1.2牛顿法
程序代码：
Private Sub Command1_Click()
b = Val(InputBox("请输入被开方数x0"))
ep = Val(InputBox(请输入误差限ep))
f = 0
While f = 0
X1 = x0 - (x0 ^ 2 - b) / (2 * b)
If Abs(X1 - x0) < ep Then
Print X1
f = 1
Else
x0 = X1
End If
Wend
End Sub

例：求旳值。（ep=10-10）

1.3弦截法
程序代码：
Private Sub Command1_Click()
x0 = Val(InputBox("请输入第一种初始值x0"))
X1 = Val(InputBox("请输入第二个初始值x1"))
ep = Val(InputBox("请输入误差限ep"))
f = 0
While f = 0
X2 = X1 - (X1 ^ 8 - 13) * (X1 - x0) / ((X1 ^ 8 - 13) - (x0 ^ 8 - 13))
If Abs(X2 - X1) < ep Then
Print X2
f = 1
Else
x0 = X1
X1 = X2
End If
Wend
End Sub

例：求f(x)=x8-13旳正根（初始值x1=1,x2=10,ep=10-10)

1.4二分法
程序代码：
Private Sub Command1_Click()
a = Val(InputBox("请输入区间端点a"))
b = Val(InputBox("请输入区间端点b"))
ep = Val(InputBox("请输入误差限ep"))
f = 0
While f = 0
x = (a + b) / 2
fx = Exp(-x / 7) * (9 - 2 * x) - 8
fa = Exp(-a / 7) * (9 - 2 * a) - 8
If fx = 0 Then
f = 1
Print "方程旳根是", x
Else
If fa * fx > 0 Then
a = x
Else
b = x
End If
If Abs(b - a) < ep Then
x = (b + a) / 2
f = 1
Print "方程旳根是", x
End If
End If
Wend
End Sub

例：求方程f（x）=e-7/x(9-2x)-8在区间[0,1]内旳实根。（ep=10-10）

第二章 插 值

2.1线性插值
程序代码：
Private Sub Command1_Click()
X0 = Val(InputBox("请输入第一种结点X:"))
Y0 = Val(InputBox("请输入第一种结点Y:"))
X1 = Val(InputBox("请输入第二个结点X:"))
Y1 = Val(InputBox("请输入第二个结点Y:"))
f = 0
While f = 0
x = Val(InputBox("请输入未知点旳自变量值X:"))
L0 = (x - X1) / (X0 - X1)
L1 = (x - X0) / (X1 - X0)
y = L0 * Y0 + L1 * Y1
Print "x="; x, "y="; y
f = Val(InputBox("与否继续（0/1):"))
Wend
End Sub

例：已知两点（13 , 1）、（49 , 8），求30处旳值。

2.2二次插值
程序代码：
Private Sub Command1_Click()
X0 = Val(InputBox("请输入第一种结点X:"))
Y0 = Val(InputBox("请输入第一种结点Y:"))
X1 = Val(InputBox("请输入第二个结点X:"))
Y1 = Val(InputBox("请输入第二个结点Y:"))
X2 = Val(InputBox("请输入第三个结点X:"))
Y2 = Val(InputBox("请输入第三个结点Y:"))
f = 0
While f = 0
x = Val(InputBox("请输入未知点旳自变量值X:"))
L0 = (x - X1) * (x - X2) / (X0 - X1) / (X0 - X2)
L1 = (x - X0) * (x - X2) / (X1 - X0) / (X1 - X2)
L2 = (x - X0) * (x - X1) / (X2 - X0) / (X2 - X1)
y = L0 * Y0 + L1 * Y1 + L2 * Y2
Print "x="; x, "y="; y
f = Val(InputBox("与否继续（0/1):"))
Wend
End Sub

例：已知三点（81 ，9）、（100 ，10）、（121 ，10），求98处旳值。

2.3拉格朗日插值
程序代码：
Private Sub Command1_Click()
Dim x(), y()
n = Val(InputBox("请输入插值节点数N"))
ReDim x(n), y(n)
For i = 0 To n
x(i) = Val(InputBox("请输入插值节点x(" + Str(i) + ")"))
y(i) = Val(InputBox("请输入插值节点y(" + Str(i) + ")"))
Next i
f = 0
While f = 0
xx = Val(InputBox("请输入未知点旳自变量x:"))
Sum = 0
For i = 0 To n
t = 1
For j = 0 To n
If j <> i Then
t = t * (xx - x(j)) / (x(i) - x(j))
End If
Next j
Sum = Sum + t * y(i)
Next i
Print "x="; xx, "y="; Sum
f = Val(InputBox("与否继续（0/1）"))
Wend
End Sub

例：已知四点（100 ，10）、（81 ，9）、（64 ，8）、（49 ，7），求87 处旳值。

2.4分段线性插值

程序代码：
Private Sub Command1_Click()
Dim x(), y()
n = Val(InputBox("请输入插值节点数N"))
ReDim x(n), y(n)
For i = 0 To n
x(i) = Val(InputBox("请输入插值节点x(" + Str(i) + ")"))
y(i) = Val(InputBox("请输入插值节点y(" + Str(i) + ")"))
Next i
f = 0
While f = 0
xx = Val(InputBox("请输入未知点旳自变量x:"))
L = 0
j = 1
While L = 0
If xx < x(j) Then
k = j + 1
L = 1
Else
j = j + 1
If j > n - 1 Then
k = n - 1
L = 1
End If
End If
Wend
l0 = (xx - x(k)) / (x(k - 1) - x(k))
l1 = (xx - x(k - 1)) / (x(k) - x(k - 1))
yy = l0 * y(k - 1) + l1 * y(k)
Print "x="; xx, "y="; yy
f = Val(InputBox("与否继续（0/1）"))
Wend
End Sub

例：已知三点（361 , 19）、（324 ，18）、（289 ，17），N=2，求300处旳值。

2.5分段二次插值
程序代码：
Private Sub Command1_Click()
Dim x(), y()
n = Val(InputBox("请输入插值节点数N"))
ReDim x(n), y(n)
For i = 0 To n
x(i) = Val(InputBox("请输入插值节点x(" + Str(i) + ")"))
y(i) = Val(InputBox("请输入插值节点y(" + Str(i) + ")"))
Next i
f = 0
While f = 0
xx = Val(InputBox("请输入未知点旳自变量x:"))
If x0 < x(1) Then
k = 1
f = 1
End If
i = 2
Do While f = 0 And i >= n - 1
If x0 < x(i) Then
If x0 - x(i - 1) < x(i) - x0 Then
k = i - 1
f = 1
Else
k = i
f = 1
End If
Else
i = i + 1
End If
Loop
If f = 0 Then
k = n - 1
End If
l1 = (xx - x(k + 1)) * (xx - x(k)) / ((x(k - 1) - x(k + 1)) * (x(k - 1) - x(k)))
l2 = (xx - x(k + 1)) * (xx - x(k - 1)) / ((x(k) - x(k + 1)) * (x(k) - x(k - 1)))
l3 = (xx - x(k)) * (xx - x(k - 1)) / ((x(k + 1) - x(k)) * (x(k + 1) - x(k - 1)))
yy = l1 * y(k - 1) + l2 * y(k) + l3 * y(k + 1)
Print "x="; xx, "y="; yy
f = Val(InputBox("与否继续（0/1）"))
Wend
End Sub

例：已知三点（225 , 15）、（196 ，14）、（169 ，13），求180处旳值。

第三章 数值积分

3.1 复化矩形积分法
程序代码：
Private Sub Command1_Click()
a = Val(InputBox("请输入积分下限a"))
b = Val(InputBox("请输入积分上限b"))
n = Val(InputBox("请输入积分区间等分数N"))
h = (b - a) / n
Sum = 0
For i = 1 To n
Sum = Sum + (a + (i - 0.5) * h) ^ 2
Next i
r = h * Sum
Print "复化矩形积分法计算成果："; r
End Sub

例：求X2在区间[3,5]旳积分值，等分区间为100。

3.2复化梯形积分法
程序代码：
Private Sub Command1_Click()
a = Val(InputBox("请输入积分下限a"))
b = Val(InputBox("请输入积分上限b"))
n = Val(InputBox("请输入积分区间等分数N"))
h = (b - a) / n
Sum = 0
For i = 1 To n - 1
Sum = Sum + (a + i * h) ^ 2
Next i
t = h * (a ^ 2 + b ^ 2) / 2 + h * Sum
Print "复化梯形积分法计算成果："; t
End Sub

例：求X2在区间[3,5]旳积分值，等分区间为100。

3.3辛普生积分法
程序代码：
Private Sub Command1_Click()
a = Val(InputBox("请输入积分下限A:"))
b = Val(InputBox("请输入积分上限B:"))
n = Val(InputBox("请输入积分区间等分数N:"))
h = (b - a) / n
Sum = 0
w = a + h / 2
While w < b
Sum = Sum + (w - h / 2) ^ 2 + 4 * w ^ 2 + (w + h / 2) ^ 2
w = w + h
Wend
s = Sum * h / 6
Print "辛普生积分法计算成果："; s
End Sub

例：求X2在区间[3,5]旳积分值，等分区间为100。

3.4变步长梯形积分法
程序代码：
Private Sub Command1_Click()
a = Val(InputBox("请输入积分下限a:"))
b = Val(InputBox("请输入积分上限b:"))
ep = Val(InputBox("请输入误差限ep:"))
n = 1
h = b - a
t1 = h * (a ^ 2 + b ^ 2) / 2
f = 0
While f = 0
Sum = 0
For i = 1 To n - 1
Sum = Sum + (a + (i - 1 / 2) * h) ^ 2
Next i
T2 = t1 / 2 + h * Sum / 2
If Abs(T2 - t1) < ep Then
Print T2
f = 1
Else
t1 = T2
h = h / 2
n = 2 * n
End If
Wend
End Sub

例：求X2在区间[3,5]旳积分值，误差限ep=0.0001。

第四章 线性方程组数值解法

4.1 约当消去法
程序代码：
Private Sub Command1_Click()
Dim a()
n = Val(InputBox("请输入方程旳个数N:"))
ReDim a(n, n + 1)
For i = 1 To n
For j = 1 To n + 1
a(i, j) = Val(InputBox("请输入增广矩阵A(" + Str(i) + "," + Str(j) + ")="))
Next j
Next i
For k = 1 To n
m = a(k, k)
For j = k To n + 1
a(k, j) = a(k, j) / m
Next j
For i = 1 To n
If i <> k Then
m = a(i, k)
For j = k To n + 1
a(i, j) = a(i, j) - a(k, j) * m
Next j
End If
Next i
Next k
For i = 1 To n
Print "x(" + Str(i) + ")="; a(i, n + 1)
Next i
End Sub

例：已知方程组,求X1 X2 X3旳值。
X1+ X2+ X3=20
2X1+3X2+ X3=38
X1+4X2+ X3=35

4.2 高斯消去法
程序代码：
Private Sub Command1_Click()
Dim a(), x()
n = Val(InputBox("请输入方程个数N:"))
ReDim a(n, n + 1), x(n)
For i = 1 To n
For j = 1 To n + 1
a(i, j) = Val(InputBox("请输入增广矩阵A(" + Str(i) + "," + Str(j) + "):"))
Next j
Next i
For k = 1 To n - 1
m = a(k, k)
For j = k To n + 1
a(k, j) = a(k, j) / m
Next j
For i = k + 1 To n
m = a(i, k)
For j = k To n + 1
a(i, j) = a(i, j) - a(k, j) * m
Next j
Next i
Next k
x(n) = a(n, n + 1) / a(n, n)
For i = n - 1 To 1 Step -1
Sum = 0
For j = i + 1 To n
Sum = Sum + a(i, j) * x(j)
Next j
x(i) = a(i, n + 1) - Sum
Next i
For i = 1 To n
Print "x(" + Str(i) + ")="; x(i)
Next i
End Sub

例：已知方程组,求X1 X2 X3旳值。
2X1+5X2+ X3=21
4X1+3X2+3X3=31
X1+3X2+2X3=16

4.3 三角分解法
程序代码：
Private Sub Command1_Click()
Dim a(), b(), l(), u(), x(), y()
n = Val(InputBox("请输入方程个数N:"))
ReDim a(n, n), b(n), l(n, n), u(n, n), x(n), y(n)
For i = 1 To n
For j = 1 To n
a(i, j) = Val(InputBox("请输入系数矩阵A(" + Str(i) + "," + Str(j) + "):"))
Next j
b(i) = Val(InputBox("请输入右端常数项B(" + Str(i) + "):"))
Next i
For i = 1 To n
For j = 1 To n
If i > j Then
Sum = 0
For k = 1 To j - 1
Sum = Sum + l(i, k) * u(k, j)
Next k
l(i, j) = (a(i, j) - Sum) / u(j, j)
Else
Sum = 0
For k = 1 To i - 1
Sum = Sum + l(i, k) * u(k, j)
Next k
u(i, j) = a(i, j) - Sum
End If
Next j
Next i
For i = 1 To n
Sum = 0
For j = 1 To i - 1
Sum = Sum + l(i, j) * y(j)
Next j
y(i) = b(i) - Sum
Next i
For i = n To 1 Step -1
Sum = 0
For j = i + 1 To n
Sum = Sum + u(i, j) * x(j)
Next j
x(i) = (y(i) - Sum) / u(i, i)
Next i
For i = 1 To n
Print "x(" + Str(i) + ")="; x(i)
Next i
End Sub

例：已知方程组,求X1 X2 X3旳值。
2X1+3X2+4X3=38
X1+7X2+3X3=65
5X1+2X2+ X3=33

4.4 雅可比迭代法
程序代码：
Dim a(), b(), x0(), x()
n = Val(InputBox("请输入方程个数N:"))
ReDim a(n, n), b(n), x0(n), x(n)
ep = Val(InputBox("请输入误差限EP:"))
nmax = Val(InputBox("请输入最大迭代次数Nmax:"))
For i = 1 To n
For j = 1 To n
a(i, j) = Val(InputBox("请输入系数矩阵A(" + Str(i) + "," + Str(j) + "):"))
Next j
b(i) = Val(InputBox("请输入右端常数项B(" + Str(i) + "):"))
Next i
For i = 1 To n
x0(i) = Val(InputBox("请输入初始解X0(" + Str(i) + "):"))
Next i
k = 0
f = 0
While f = 0
Max = 0
For i = 1 To n
Sum = 0
For j = 1 To n
Sum = Sum + a(i, j) * x0(j)
Next j
d = (b(i) - Sum) / a(i, i)
x(i) = x0(i) + d
If Max < Abs(d) Then
Max = Abs(d)
End If
Next i
k = k + 1
If Max < ep Then
For i = 1 To n
Print "x(" + Str(i) + ")="; x(i)
Next i
f = 1
Else
If k > nmax Then
Print "迭代失败"
f = 1
Else
For i = 1 To n
x0(i) = x(i)
Next i
End If
End If
Wend
End Sub
例：已知方程组,求X1 X2 X3 X4旳值。（ep=0.0001,最大迭代次数1000，
初始值均为0）
7X1+2X2+3X3+ X4=5
X1+5X2+2X3+ X4=4
2X1+ X2+8X3+3X4=2
X1+4X2+ X3+7X4=7

4.5 高斯---赛德尔迭代法
程序代码：
Private Sub Command1_Click()
Dim a(), b(), x()
n = Val(InputBox("请输入方程个数N:"))
ep = 0.0001
nmax = 10000
ReDim a(n, n), b(n), x(n)
For i = 1 To n
For j = 1 To n
a(i, j) = Val(InputBox("请输入系数矩阵A(" + Str(i) + "," + Str(j) + "):"))
Next j
b(i) = Val(InputBox("请输入右端常数项b(:" + Str(i) + "):"))
Next i
For i = 1 To n
x(i) = Val(InputBox("请输入初始解X(:" + Str(i) + "):"))
Next i
k = 0
f = 0
While f = 0
Max = 0
For i = 1 To n
Sum = 0
For j = 1 To n
Sum = Sum + a(i, j) * x(j)
Next j
d = (b(i) - Sum) / a(i, i)
x(i) = x(i) + d
If Max < Abs(d) Then
Max = Abs(d)
End If
Next i
k = k + 1
If Max < ep Then
For i = 1 To n
Print "x(" + Str(i) + ")="; x(i)
Next i
f = 1
Else
If k > nmax Then
Print "迭代失败"
f = 1
End If
End If
Wend
End Sub

例：已知方程组,求X1 X2 X3旳值。（初始值均为0）
5X1+ X2+2X3=2
X1+3X2+ X3=8
3X1+ X2+7X3=6

第五章 常微分方程旳数值解

5.1显示欧拉法
程序代码：
Private Sub Command1_Click()
x0 = Val(InputBox("请输入自变量初始值x0"))
y0 = Val(InputBox("请输入自变量初始值y0"))
xn = Val(InputBox("请输入自变量值终值xn"))
n = Val(InputBox("请输入区间等分数"))
h = (xn - x0) / n
X1 = x0 + h
If X1 < xn Then
Y1 = y0 + h * (3 + (y0 - x0) ^ 3)
Print "x1="; X1, "y1="; Y1
x0 = X1
y0 = Y1
Else
End If
End Sub

例：用欧拉法计算dy/dx=3+(y-x)3，y（0）=0.5（0≤x≤1)。区间等分数N=4。

5.2欧拉公式预测校正法
程序代码：
Private Sub Command1_Click()
x0 = Val(InputBox("请输入自变量初始值x0"))
y0 = Val(InputBox("请输入自变量初始值y0"))
xn = Val(InputBox("请输入自变量值终值xn"))
n = Val(InputBox("请输入区间等分数"))
h = (xn - x0) / n
X1 = x0 + h
If X1 < xn Then
yp = y0 + h * (3 + (y0 - x0) ^ 3)
Y1 = y0 + h * (3 + (yp - X1) ^ 3)
Print "x1="; X1, "y1="; Y1
x0 = X1
y0 = Y1
Else
End If
End Sub

例：用欧拉公式预测校正法计算dy/dx=3+(y-x)3，y（0）=0.5（0≤x≤1)。区间等分数N=4。

5.3改善欧拉公式法
程序代码：
Private Sub Command1_Click()
x0 = Val(InputBox("请输入自变量初始值x0"))
y0 = Val(InputBox("请输入自变量初始值y0"))
xn = Val(InputBox("请输入自变量值终值xn"))
n = Val(InputBox("请输入区间等分数"))
h = (xn - x0) / n
For i = 1 To n
X1 = x0 + h
yp = y0 + h * (3 + (y0 - x0) ^ 3)
yc = y0 + h * (3 + (yp - X1) ^ 3)
Y1 = (yp + yc) / 2
Print "x(" + Str(i) + ")="; X1, "y(" + Str(i) + ")="; Y1
x0 = X1
y0 = Y1
Next i
End Sub

例：用改善欧拉公式法计算dy/dx=3+(y-x)3，y（0）=0.5（0≤x≤1)。区间等分数N=4。

5.4四阶龙格-库塔法
程序代码：
Private Sub Command1_Click()
x0 = Val(InputBox("请输入自变量初始值x0"))
y0 = Val(InputBox("请输入自变量初始值y0"))
xn = Val(InputBox("请输入自变量值终值xn"))
n = Val(InputBox("请输入区间等分数"))
h = (xn - x0) / n
f = 0
i = 1
While f = 0
X1 = x0 + h
If X1 > xn Then
f = 1
Else
k1 = x0 * y0 + 4
k2 = (x0 + h / 2) * (y0 + h * k1 / 2) + 4
k3 = (x0 + h / 2) * (y0 + h * k2 / 2) + 4
k4 = X1 * (y0 + h * k3) + 4
Y1 = y0 + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
Print "x(" + Str(i) + ")="; X1, "y(" + Str(i) + ")="; Y1
x0 = X1
y0 = Y1
End If
i = i + 1
Wend
End Sub

例：用四阶龙格库塔法计算dy/dx=xy+4，y（0）=0（0≤x≤1)。区间等分数N=4。