Skip to content

Instantly share code, notes, and snippets.

@baishuai

baishuai/NA数值分析实验.md

Last active June 14, 2016 23:19
Show Gist options
  • Save baishuai/fd88a0fe0a4301081c27 to your computer and use it in GitHub Desktop.
Save baishuai/fd88a0fe0a4301081c27 to your computer and use it in GitHub Desktop.
Download ZIP
数值分析实验2014年 蔡懿慈
Raw
lab1.go
/*
* @Author: Bai Shuai
* @Date: 2014-05-12 13:53:54
* @Last Modified by: bai
* @Last Modified time: 2014-05-22 22:49:14
*/
/* 运行结果
截断误差与舍入误差 "精确ln2值:" 0.693147182465
误差限制:5.000000e-06 n=40871 ln2=0.693152129650 dx=4.9471855e-06
误差限制:5.000000e-07 n=61343 ln2=0.693147659302 dx=4.7683716e-07
*/
package main
import "fmt"
const (
ln2value float32 = 0.693147180560
epsilon1 float32 = 0.5e-5 //1/2 x 10^-5
epsilon2 float32 = 0.5e-6 //1/1 x 10^-6
)
//abs 绝对值函数
func abs(x float32) float32 {
switch {
case x < 0:
return -x
case x == 0:
return 0 // return correctly abs(-0)
}
return x
}
//ln2 根据级数公式求ln2值，根据级数n-1时的ln2值求级数n时的ln2值
func ln2(n int, value float32) float32 {
if n&1 == 1 {
value += 1.0 / float32(n)
} else {
value -= 1.0 / float32(n)
}
return value
}
//ln2Accu 计算满足给定精度的ln2值，返回n值与计算的ln2值
//para：
// n 级数公式的n值
// value 级数公式在n-1时求到的ln2值
// accu 要求的精度
func ln2Accu(n int, value, accu float32) (int, float32, float32) {
value = ln2(n, value)
for abs(value-ln2value) >= accu {
n++
value = ln2(n, value)
}
return n, value, abs(value - ln2value)
}
func main() {
fmt.Printf("截断误差与舍入误差 \"精确ln2值:\" %16.12f\n", ln2value)
n, v, d := ln2Accu(1, 0.0, epsilon1)
fmt.Printf("误差限制:%e n=%d ln2=%.12f dx=%v\n", epsilon1, n, v, d)
n, v, d = ln2Accu(n+1, v, epsilon2)
fmt.Printf("误差限制:%e n=%d ln2=%.12f dx=%v\n", epsilon2, n, v, d)
}
Raw
lab2.go
/*
* @Author: Bai Shuai
* @Date: 2014-05-15 10:49:42
* @Last Modified by: Bai Shuai
* @Last Modified time: 2014-05-15 11:01:13
*/
package main
import "fmt"
func lab2f1(x float64) float64 {
return 1.0 / (1 + 16*x*x)
}
func main() {
}
Raw
lab3.go
/*
* @Author: Bai Shuai
* @Date: 2014-05-11 20:56:50
* @Last Modified by: bai
* @Last Modified time: 2014-05-31 22:45:41
*/
/* 运行结果
-27.121212122148506+48.88720538728549X^1-0.5111111111132265X^2+0.0071380471380646324X^3
20 805 803.2828282827878 1.7171717172121816
25 985 987.1464646464822 2.1464646464821726
30 1170 1172.2222222222574 2.2222222222574146
35 1365 1363.8636363636624 1.1363636363375917
40 1570 1567.424242424245 2.5757575757550057
45 1790 1788.2575757575548 1.7424242424451677
50 2030 2031.717171717139 1.717171717139081
55 2300 2303.156565656547 3.1565656565471727
60 2610 2607.9292929293265 2.0707070706735067
*/
package main
import (
"fmt"
"math"
"strconv"
)
var (
mSize, nSize int //m为点集数，n为拟合曲线次数
xVal, omega []float64 //拟合数据以及权重
fx map[float64]float64
aStar []float64 //拟合曲线系数
alpha, beta []float64 //计算Pk需要的中间值
pFunc []*polyFunc
)
type polyFunc struct {
rank int
factor []float64
}
func (p *polyFunc) String() (str string) {
str = strconv.FormatFloat(p.factor[0], 'f', -1, 64)
for i := 1; i < p.rank; i++ {
if p.factor[i] > 0 {
str += "+"
}
str += strconv.FormatFloat(p.factor[i], 'f', -1, 64)
str += "X^" + strconv.Itoa(i)
}
return str
}
func (p *polyFunc) Val() func(float64) float64 {
return func(x float64) (val float64) {
for i := 0; i < p.rank; i++ {
val += p.factor[i] * math.Pow(x, float64(i))
}
return val
}
}
func (p *polyFunc) Mult(q *polyFunc) (s *polyFunc) {
s = new(polyFunc)
s.rank = p.rank + q.rank - 1
s.factor = make([]float64, s.rank)
for i := 0; i < p.rank; i++ {
for j := 0; j < q.rank; j++ {
s.factor[i+j] += p.factor[i] * q.factor[j]
}
}
return s
}
func (p *polyFunc) Add(q *polyFunc) (s *polyFunc) {
s = new(polyFunc)
if p.rank > q.rank {
p, q = q, p
}
s.rank = q.rank
s.factor = make([]float64, s.rank)
for i := 0; i < p.rank; i++ {
s.factor[i] = p.factor[i] + q.factor[i]
}
for i := p.rank; i < q.rank; i++ {
s.factor[i] = q.factor[i]
}
if s.factor[s.rank-1] == 0 {
s.rank--
}
return s
}
func (p *polyFunc) Neg() *polyFunc {
for i := 0; i < p.rank; i++ {
p.factor[i] = -p.factor[i]
}
return p
}
func innerProduct(fa, fb func(float64) float64) (ip float64) {
for i := 0; i < mSize; i++ {
ip += fa(xVal[i]) * fb(xVal[i]) * omega[i]
}
return ip
}
func Solve(x, y, om []float64, n int) (func(float64) float64, string) {
xVal = x
mSize = len(xVal)
nSize = n
if om == nil {
omega = make([]float64, mSize)
for i := 0; i < mSize; i++ {
omega[i] = 1.0
}
} else {
omega = om
}
fx = make(map[float64]float64)
for i := 0; i < mSize; i++ {
fx[xVal[i]] = y[i]
}
alpha = make([]float64, nSize)
beta = make([]float64, nSize)
pFunc = make([]*polyFunc, nSize)
aStar = make([]float64, nSize)
for i := 0; i < nSize; i++ {
pFunc[i] = calcPunck(i)
aStar[i] = innerProduct(func(x float64) float64 {
return fx[x]
}, pFunc[i].Val()) / innerProduct(pFunc[i].Val(), pFunc[i].Val())
}
fans := &polyFunc{1, []float64{0}}
for i := 0; i < nSize; i++ {
fans = pFunc[i].Mult(&polyFunc{1, []float64{aStar[i]}}).Add(fans)
}
return fans.Val(), fans.String()
}
func calcAlpha(k int) float64 {
return innerProduct(func(x float64) float64 {
return x * pFunc[k-1].Val()(x)
}, pFunc[k-1].Val()) / innerProduct(pFunc[k-1].Val(), pFunc[k-1].Val())
}
func calcBeta(k int) float64 {
return innerProduct(pFunc[k].Val(), pFunc[k].Val()) /
innerProduct(pFunc[k-1].Val(), pFunc[k-1].Val())
}
func calcPunck(i int) (p *polyFunc) {
if i == 0 {
p = &polyFunc{1, []float64{1.0}}
} else if i == 1 {
alpha[i] = calcAlpha(i)
p = pFunc[0].Mult(&polyFunc{2, []float64{-alpha[i], 1}})
} else {
alpha[i] = calcAlpha(i)
beta[i-1] = calcBeta(i - 1)
p = pFunc[i-1].Mult(&polyFunc{2, []float64{-alpha[i], 1}}).Add(
pFunc[i-2].Mult(&polyFunc{1, []float64{-beta[i-1]}}))
}
return p
}
func main() {
x := []float64{20, 25, 30, 35, 40, 45, 50, 55, 60}
y := []float64{805, 985, 1170, 1365, 1570, 1790, 2030, 2300, 2610}
fAns, fstr := Solve(x, y, nil, 4)
fmt.Println(fstr)
for i := 0; i < len(x); i++ {
fmt.Println(x[i], y[i], fAns(x[i]), math.Abs(fAns(x[i])-y[i]))
}
}
Raw
lab4.go
/*
* @Author: Bai Shuai
* @Date: 2014-05-14 21:12:12
* @Last Modified by: bai
* @Last Modified time: 2014-06-02 16:53:49
*/
/* 运行结果
Ln(2)
复合辛普森 0.6931471820207081 1.4607628440188591e-09
[0.75]
[0.7083333333333333 0.6944444444444443]
[0.6970238095238095 0.6932539682539682 0.6931746031746032]
[0.6941218503718504 0.6931545306545307 0.6931479014812348 0.6931474776448321]
[0.6933912022075268 0.6931476528194189 0.6931471942970782 0.6931471830719329 0.6931471819167451]
[0.693208208269249 0.6931472102898231 0.69314718078785 0.6931471805734177 0.6931471805636197 0.693147180562297]
Romberg外推 0.693147180562297 2.3516744107610066e-12
复合高斯 0.6931471793189892 -1.2409561156445648e-09
Pi
复合辛普森 3.141592653552836 -3.695710404372221e-11
[3]
[3.1 3.1333333333333337]
[3.131176470588235 3.1415686274509804 3.1421176470588232]
[3.138988494491089 3.141592502458707 3.141594094125889 3.141585783761874]
[3.140941612041389 3.1415926512248227 3.1415926611425635 3.1415926383967965 3.141592665277718]
[3.1414298931749745 3.141592653552836 3.1415926537080368 3.1415926535900285 3.1415926536496097 3.1415926536382432]
[3.141551963485656 3.1415926535892167 3.1415926535916423 3.141592653589795 3.141592653589794 3.1415926535897354 3.1415926535897234]
Romberg外推 3.1415926535897234 -6.972200594645983e-14
复合高斯 3.141592653681139 9.134604184168893e-11
*/
package main
import (
"fmt"
"math"
)
func fLn2(x float64) float64 {
return 1.0 / x
}
func fPi(x float64) float64 {
return 4.0 / (x*x + 1)
}
// ComSimpson 复合辛普森求积公式
func ComSimpson(f func(float64) float64, a, b float64, n int) float64 {
h := (b - a) / float64(n)
ans := f(a) + f(b) + 4*f(a+h/2)
var xk float64
for i := 1; i < n; i++ {
xk = a + h*float64(i)
ans = ans + 4*f(xk+h/2) + 2*f(xk)
}
ans = ans * h / 6
return ans
}
// Romberg 外推求积分
func Romberg(f func(float64) float64, a, b, accu float64) float64 {
h := b - a
T := make([][]float64, 0, 5)
k := 0
T = append(T, []float64{(f(a) + f(b)) * h / 2})
fmt.Println(T[0])
for {
tmp := 0.0
for i := 0; i < 1<<uint(k); i++ {
tmp += f(a + h*float64(i) + h/2)
}
k++
h = h / 2
Tk := make([]float64, k+1)
Tk[0] = h*tmp + T[k-1][0]/2
for i := 1; i <= k; i++ {
tmp = math.Pow(4.0, float64(i))
Tk[i] = Tk[i-1]*tmp/(tmp-1.0) - T[k-1][i-1]/(tmp-1.0)
}
T = append(T, Tk)
fmt.Println(T[k])
if math.Abs(T[k][k]-T[k-1][k-1]) < accu {
break
}
}
return T[k][k]
}
// ComGauss 求积分
func ComGauss(f func(float64) float64, a, b float64, n int) float64 {
h := (b - a) / float64(n)
h2sq3 := h / (2.0 * math.Sqrt(3.0))
var xk, ans float64
for i := 0; i < n; i++ {
xk = a + h*float64(i) + h/2.0
ans += f(xk+h2sq3) + f(xk-h2sq3)
}
return ans * h / 2.0
}
func main() {
fmt.Println("Ln(2)")
ln2 := ComSimpson(fLn2, 1, 2, 34)
fmt.Println("复合辛普森", ln2, ln2-math.Ln2)
ln2 = Romberg(fLn2, 1, 2, 0.5e-8)
fmt.Println("Romberg外推", ln2, ln2-math.Ln2)
ln2 = ComGauss(fLn2, 1, 2, 32)
fmt.Println("复合高斯", ln2, ln2-math.Ln2)
fmt.Println("\nPi")
pi := ComSimpson(fPi, 0, 1, 16)
fmt.Println("复合辛普森", pi, pi-math.Pi)
pi = Romberg(fPi, 0, 1, 0.5e-8)
fmt.Println("Romberg外推", pi, pi-math.Pi)
pi = ComGauss(fPi, 0, 1, 13)
fmt.Println("复合高斯", pi, pi-math.Pi)
}
Raw
lab5.go
/*
* @Author: Bai Shuai
* @Date: 2014-05-11 20:42:29
* @Last Modified by: bai
* @Last Modified time: 2014-05-31 22:47:42
*/
//非线性方程求根
//二分法、牛顿法、迭代法
//1. 求解非线性方程
//2. 研究迭代函数、迭代初值对函数收敛性及收敛速度的影响
/* 实验结果
迭代次数 1, 迭代值 0.0000000000
迭代次数 2, 迭代值 1.5000000000
迭代次数 3, 迭代值 0.7500000000
迭代次数 4, 迭代值 0.3750000000
迭代次数 5, 迭代值 0.1875000000
迭代次数 6, 迭代值 0.2812500000
迭代次数 7, 迭代值 0.3281250000
迭代次数 8, 迭代值 0.3515625000
迭代次数 9, 迭代值 0.3398437500
迭代次数 10, 迭代值 0.3457031250
迭代次数 11, 迭代值 0.3427734375
迭代次数 12, 迭代值 0.3442382812
迭代次数 13, 迭代值 0.3449707031
迭代次数 14, 迭代值 0.3453369141
迭代次数 15, 迭代值 0.3455200195
迭代次数 16, 迭代值 0.3456115723
迭代次数 17, 迭代值 0.3456573486
迭代次数 18, 迭代值 0.3456344604
迭代次数 19, 迭代值 0.3456230164
迭代次数 20, 迭代值 0.3456287384
迭代次数 21, 迭代值 0.3456258774
迭代次数 22, 迭代值 0.3456273079
迭代次数 23, 迭代值 0.3456280231
迭代次数 24, 迭代值 0.3456276655
-------------------我是分割线----------------
Newton迭代法, 初值0.0000000000
迭代次数 1, 迭代值 0.3333333333
迭代次数 2, 迭代值 0.3456790123
迭代次数 3, 迭代值 0.3456273932
-------------------我是分割线----------------
迭代法 初值 0.0000000000, 次数 10
迭代次数 1, 迭代值 0.7937005260
迭代次数 2, 迭代值 0.9643617579
迭代次数 3, 迭代值 0.9940246594
迭代次数 4, 迭代值 0.9990031165
迭代次数 5, 迭代值 0.9998338251
迭代次数 6, 迭代值 0.9999723034
迭代次数 7, 迭代值 0.9999953839
迭代次数 8, 迭代值 0.9999992306
迭代次数 9, 迭代值 0.9999998718
迭代次数 10, 迭代值 0.9999999786
-------------------我是分割线----------------
迭代法 初值 0.0000000000, 次数 10
迭代次数 1, 迭代值 -1.0000000000
迭代次数 2, 迭代值 -3.0000000000
迭代次数 3, 迭代值 -55.0000000000
迭代次数 4, 迭代值 -332751.0000000000
迭代次数 5, 迭代值 -73686529681121504.0000000000
迭代次数 6, 迭代值 -800192186653981478358652672830638575390678549790720.0000000000
迭代次数 7, 迭代值 -1024738174056893932527762049050938485592367730189188160409266654928227462101562808957796957393014967922717761885876347463950283744054040677931592393424896.0000000000
迭代次数 8, 迭代值 -Inf
迭代次数 9, 迭代值 -Inf
迭代次数 10, 迭代值 -Inf
-------------------我是分割线----------------
Newton迭代法, 初值1.5000000000
迭代次数 1, 迭代值 1.3478260870
迭代次数 2, 迭代值 1.3252003990
迭代次数 3, 迭代值 1.3247181740
-------------------我是分割线----------------
Newton迭代法, 初值0.0000000000
迭代次数 1, 迭代值 -1.0000000000
迭代次数 2, 迭代值 -0.5000000000
迭代次数 3, 迭代值 -3.0000000000
迭代次数 4, 迭代值 -2.0384615385
迭代次数 5, 迭代值 -1.3902821472
迭代次数 6, 迭代值 -0.9116118977
迭代次数 7, 迭代值 -0.3450284967
迭代次数 8, 迭代值 -1.4277507040
迭代次数 9, 迭代值 -0.9424179125
迭代次数 10, 迭代值 -0.4049493572
迭代次数 11, 迭代值 -1.7069046452
迭代次数 12, 迭代值 -1.1557563611
迭代次数 13, 迭代值 -0.6941918133
迭代次数 14, 迭代值 0.7424942987
迭代次数 15, 迭代值 2.7812959407
迭代次数 16, 迭代值 1.9827252470
迭代次数 17, 迭代值 1.5369273798
迭代次数 18, 迭代值 1.3572624832
迭代次数 19, 迭代值 1.3256630944
迭代次数 20, 迭代值 1.3247187886
迭代次数 21, 迭代值 1.3247179572
*/
package main
import (
"fmt"
"math"
)
const (
accuracy = 1e-6
)
func f1(x float64) float64 {
return math.Pow(x, 3)*2 - math.Pow(x, 2) + 3*x - 1
}
func df1(x float64) float64 {
return math.Pow(x, 2)*6 - 2*x + 3
}
// BinCalc 二分法解非线性方程
func BinCalc(left, right float64, f func(float64) float64, accu float64) {
var mid, fa, fmid float64
var k int
fa = f(left)
// fb = f(right)
for {
mid = (left + right) / 2
fmid = f(mid)
k++
fmt.Printf("迭代次数 %d, 迭代值 %.10f\n", k, mid)
if right-left < accu {
break
} else {
if fmid*fa < 0 { //fa的值不需要更新
right = mid
} else {
left = mid
}
}
}
}
func iterf1(x float64) float64 {
return math.Pow((x+1)/2, 1.0/3)
}
func iterf2(x float64) float64 {
return math.Pow(x, 3)*2 - 1
}
// Iteration 实现了不动点迭代法解非线性方程
func Iteration(init float64, f func(float64) float64, times int) {
fmt.Printf("迭代法 初值 %.10f, 次数 %d\n", init, times)
n := 0
for n < times {
n++
init = f(init)
fmt.Printf("迭代次数 %d, 迭代值 %.10f\n", n, init)
}
}
func f2b(x float64) float64 {
return math.Pow(x, 3) - x - 1
}
func df2b(x float64) float64 {
return math.Pow(x, 2)*3 - 1
}
// Newton 实现了牛顿法解非线性方程
func Newton(init float64, f, df func(float64) float64, accu float64) {
fmt.Printf("Newton迭代法, 初值%.10f\n", init)
k := 0
xk, xk1 := init, init
var fxk, dfxk, delta float64
fxk, dfxk = f(xk), df(xk)
for {
k++
xk = xk1
xk1 = xk - fxk/dfxk
fmt.Printf("迭代次数 %d, 迭代值 %.10f\n", k, xk1)
fxk = f(xk1)
dfxk = df(xk1)
//检查精度
if math.Abs(xk1) < 1 {
delta = math.Abs(xk1 - xk)
} else {
delta = math.Abs((xk1 - xk) / xk1)
}
if delta < accu || math.Abs(fxk) < accu {
break
}
}
}
// 打印实验信息
func init() {
fmt.Println(`实验五-非线性方程求根
二分法、牛顿法、迭代法
1. 求解非线性方程
2. 研究迭代函数、迭代初值对函数收敛性及收敛速度的影响`)
}
func main() {
BinCalc(-3.0, 3.0, f1, accuracy)
fmt.Println("-------------------我是分割线----------------")
Newton(0.0, f1, df1, accuracy)
fmt.Println("-------------------我是分割线----------------")
Iteration(0.0, iterf1, 10)
fmt.Println("-------------------我是分割线----------------")
Iteration(0.0, iterf2, 10)
fmt.Println("-------------------我是分割线----------------")
Newton(1.5, f2b, df2b, accuracy)
fmt.Println("-------------------我是分割线----------------")
Newton(0.0, f2b, df2b, accuracy)
}
Raw
lab6.go
/*
* @Author: Bai Shuai
* @Date: 2014-05-12 18:45:55
* @Last Modified by: bai
* @Last Modified time: 2014-05-31 16:12:03
*/
/* 运行结果
Hilbert(3)条件数 748
Hilbert(4)条件数 28374.999999999996
[0.9999999987548338 1.0000001067849733 0.9999977378614773 1.0000204794185086 0.9999026418473493 1.0002669070133299 0.9995630884702723 1.0004214011575985 0.9997791407962119 1.0000484987218523]
0.00040664622819364116 0.00043691152972769043
引入微小偏差
[1.0000099979202124 0.9995051783319866 1.0079162240487853 0.9399741640249581 1.252089684043749 0.36981472847701546 1.9602324057847962 0.12554140511875153 1.4372125629498806 0.9077026506617476]
1e-7 0.8937161403019374 0.9602324057847962
[0.9999999987517899 1.0000001068204887 0.9999977407626763 1.0000204270666115 0.999902990991373 1.0002657230669296 0.9995653379322409 1.000418976008124 0.9997805271605786 1.0000481722546402]
0.0004045525887202972 0.0004346620677591062
LU分解法与改进平方根分解法各运行100000次的时间
LU 488.518679ms
平方根 403.02228ms
增大希尔伯特矩阵的阶，观察误差
[0.9999999999999981 1.0000000000000122 0.9999999999999875]
3 1.3322676295501878e-14 1.2545520178264269e-14
[0.9999999999999775 1.0000000000002485 0.999999999999402 1.0000000000003886]
4 5.681011217006926e-13 5.979661210631093e-13
[0.9999999999999655 1.0000000000005238 0.9999999999980628 1.0000000000025913 0.9999999999988525]
5 2.292166456641098e-12 2.5912605394751154e-12
[0.999999999999228 1.000000000021937 0.9999999998517923 1.0000000003853693 0.999999999574584 1.00000000016768]
6 3.836870821061211e-10 4.254160357319847e-10
[0.9999999999944524 1.0000000002215987 0.9999999978643687 1.0000000083033338 0.9999999847780691 1.0000000131521882 0.9999999956820236]
7 1.4256008373791929e-08 1.5221930937947548e-08
[0.9999999999662691 1.00000000180906 0.9999999763726758 1.0000001278681026 0.9999996557641158 1.0000004870421637 0.9999996534271247 1.0000000977747472]
8 4.3677961769628126e-07 4.870421637104272e-07
[0.9999999997602121 1.0000000164521592 0.9999997226859642 1.0000019734507875 0.9999927795056808 1.0000147132093888 0.9999831308861128 1.0000101748012682 0.9999974890861828]
9 1.4646253767125472e-05 1.686911388720791e-05
[0.9999999987548338 1.0000001067849733 0.9999977378614773 1.0000204794185086 0.9999026418473493 1.0002669070133299 0.9995630884702723 1.0004214011575985 0.9997791407962119 1.0000484987218523]
10 0.00040664622819364116 0.00043691152972769043
[0.9999999947511966 1.0000005467463535 0.9999858683437002 1.0001575494686312 0.9990635370043287 1.0032863331278048 0.9928557892293705 1.0097264868815559 0.9919301559258117 1.0037298503490202 0.999263885025643]
11 0.008775731886447868 0.009726486881555862
[0.9999999780153819 1.0000027186586886 0.9999161379832394 1.0011249174234058 0.9918589744462714 1.0353847131912473 0.9023148128133737 1.1754194756694725 0.7957677562849788 1.1486625704207052 0.9385254722941213 1.011022484909254]
12 0.1797003800983541 0.20423224371502124
[1.00000010084578 0.9999841078014012 1.000614076745965 0.9897633935558622 1.0919750464302016 0.5010034238206807 2.740851874285045 -3.0360583613106233 7.283884268701919 -5.493529661418356 5.270874315267408 -0.618191404610789 1.2688289098048169]
13 5.828312474054975 6.493529661418356
*/
package main
import (
"fmt"
"math"
"time"
)
// Hilbert 返回n阶Hilbert矩阵
func Hilbert(n int) [][]float64 {
mxh := make([][]float64, n)
for i := 0; i < n; i++ {
mxh[i] = make([]float64, n)
for j := 0; j < n; j++ {
mxh[i][j] = 1.0 / float64(i+j+1)
}
}
return mxh
}
// Comb 计算组合数C(m,n)
func Comb(m, n int) int {
//=m!/((m-n)!*n!)
if n == 0 {
return 1
}
mnf := n + 1
for i := n + 2; i <= m; i++ {
mnf *= i
}
msnf := 1
for i := 2; i <= m-n; i++ {
msnf *= i
}
return mnf / msnf
}
// InHilbert 返回n阶Hilbert矩阵的逆矩阵
func InHilbert(n int) [][]float64 {
imxh := make([][]float64, n)
for i := 0; i < n; i++ {
imxh[i] = make([]float64, n)
for j := i; j < n; j++ {
imxh[i][j] = float64((i + j + 1) * Comb(n+i, n-j-1) *
Comb(n+j, n-i-1) * Comb(i+j, i) * Comb(i+j, i))
if (i+j)%2 == 1 {
imxh[i][j] = -imxh[i][j]
}
}
}
for i := 1; i < n; i++ {
for j := 0; j < i; j++ {
imxh[i][j] = imxh[j][i]
}
}
return imxh
}
// NormInf 两个向量误差的无穷范数
func NormInf(x1, x2 []float64) (xo float64) {
for i := 0; i < len(x1); i++ {
if math.Abs(x1[i]-x2[i]) > xo {
xo = math.Abs(x1[i] - x2[i])
}
}
return xo
}
// NormMInf 矩阵m的无穷范数
func NormMInf(mx [][]float64) (x float64) {
for i := 0; i < len(mx); i++ {
tmp := 0.0
for j := 0; j < len(mx); j++ {
tmp += math.Abs(mx[i][j])
}
if tmp > x {
x = tmp
}
}
return x
}
// CondInfWithInverse 计算条件数
func CondInfWithInverse(mx, imx [][]float64) float64 {
return NormMInf(mx) * NormMInf(imx)
}
// A2LL 对矩阵进行LL分解（改进平方根法）
func A2LL(mx [][]float64) [][]float64 {
n := len(mx)
for i := 1; i < n; i++ {
for j := 0; j < i; j++ {
// mx[j][i] = mx[i][j]
// tmp := mx[j][i]
for k := 0; k < j; k++ {
mx[j][i] -= mx[k][i] * mx[j][k]
}
// mx[j][i] = tmp
}
for j := 0; j < i; j++ {
mx[i][j] = mx[j][i] / mx[j][j]
}
//tmp := mx[i][i]
for k := 0; k < i; k++ {
mx[i][i] -= mx[k][i] * mx[i][k]
}
//mx[i][i] = tmp
}
return mx
}
// LLUb2x 使用平方根分解求解线性方程组 mxa*x = b
func LLUb2x(mxa [][]float64, b []float64) []float64 {
n := len(b)
ll := A2LL(mxa)
y := make([]float64, n)
y[0] = b[0]
for i := 1; i < n; i++ {
tmp := b[i]
for k := 0; k < i; k++ {
tmp -= ll[i][k] * y[k]
}
y[i] = tmp
}
x := make([]float64, n)
x[n-1] = y[n-1] / ll[n-1][n-1]
for i := n - 2; i >= 0; i-- {
tmp := y[i] / ll[i][i]
for k := i + 1; k < n; k++ {
tmp -= ll[k][i] * x[k]
}
x[i] = tmp
}
return x
}
// A2LU 对矩阵进行LU分解
func A2LU(mx [][]float64) [][]float64 {
n := len(mx)
//计算U第1行，L第1列
for i := 1; i < n; i++ {
mx[i][0] = mx[i][0] / mx[0][0]
}
for r := 1; r < n; r++ {
for i := r; i < n; i++ {
for k := 0; k < r; k++ {
mx[r][i] -= mx[r][k] * mx[k][i]
}
}
for i := r + 1; i < n; i++ {
for k := 0; k < r; k++ {
mx[i][r] -= mx[i][k] * mx[k][r]
}
mx[i][r] /= mx[r][r]
}
}
return mx
}
// LUb2x 使用LU分解法求解方程组
func LUb2x(mxa [][]float64, b []float64) []float64 {
n := len(b)
lu := A2LU(mxa)
y := make([]float64, n)
y[0] = b[0]
for i := 1; i < n; i++ {
y[i] = b[i]
for k := 0; k < i; k++ {
y[i] -= lu[i][k] * y[k]
}
}
x := make([]float64, n)
x[n-1] = y[n-1] / lu[n-1][n-1]
for i := n - 2; i >= 0; i-- {
x[i] = y[i]
for k := i + 1; k < n; k++ {
x[i] -= lu[i][k] * x[k]
}
x[i] /= lu[i][i]
}
return x
}
func martixMulLiner(hi [][]float64, x []float64) []float64 {
n := len(x)
b := make([]float64, n)
for i := 0; i < n; i++ {
for j := 0; j < n; j++ {
b[i] += hi[i][j] * x[i]
}
}
return b
}
func lab62(n int) (r, dx float64) {
hn := Hilbert(n)
x := make([]float64, n)
for i := 0; i < n; i++ {
x[i] = 1.0
}
b := martixMulLiner(hn, x)
xhat := LUb2x(hn, b)
fmt.Println(xhat)
r = NormInf(b, martixMulLiner(Hilbert(n), xhat))
dx = NormInf(x, xhat)
fmt.Println(n, r, dx)
return
}
func main() {
fmt.Println("Hilbert(3)条件数", CondInfWithInverse(Hilbert(3), InHilbert(3)))
fmt.Println("Hilbert(4)条件数", CondInfWithInverse(Hilbert(4), InHilbert(4)))
x := []float64{1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
b := martixMulLiner(Hilbert(10), x)
//h10 = A2LU(h10)
xhat := LUb2x(Hilbert(10), b)
fmt.Println(xhat)
r := NormInf(b, martixMulLiner(Hilbert(10), xhat))
dx := NormInf(x, xhat)
fmt.Println(r, dx)
b[0] = b[0] + 1e-7
xhat = LUb2x(Hilbert(10), b)
fmt.Println(xhat)
r = NormInf(b, martixMulLiner(Hilbert(10), xhat))
dx = NormInf(x, xhat)
fmt.Println("1e-7", r, dx)
fmt.Println()
x = []float64{1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
b = martixMulLiner(Hilbert(10), x)
xhat = LLUb2x(Hilbert(10), b)
fmt.Println(xhat)
r = NormInf(b, martixMulLiner(Hilbert(10), xhat))
dx = NormInf(x, xhat)
fmt.Println(r, dx)
fmt.Printf("\nLU分解法与改进平方根分解法各运行100000次的时间\n")
now := time.Now()
for i := 0; i < 100000; i++ {
LUb2x(Hilbert(10), b)
}
fmt.Println(time.Now().Sub(now).String())
now = time.Now()
for i := 0; i < 100000; i++ {
LLUb2x(Hilbert(10), b)
}
fmt.Println(time.Now().Sub(now).String())
fmt.Printf("\n增大希尔伯特矩阵的阶，观察误差\n")
for i := 3; i < 25; i++ {
r, dx := lab62(i)
if r > 1 || dx > 1 {
break
}
}
}
Raw
lab7.go
/*
* @Author: Bai Shuai
* @Date: 2014-05-12 13:57:18
* @Last Modified by: Bai Shuai
* @Last Modified time: 2014-05-14 18:09:03
*/
//Iterative Methods for Sparse Linear Systems
package main
import (
"fmt"
"math"
)
var (
n = 100
h = 1.0 / 100
)
func matrixA(epsilon float64) [][]float64 {
mA := make([][]float64, n)
eh := epsilon + h
n2eh := -epsilon*2.0 - h
mA[0] = make([]float64, n)
mA[0][0] = n2eh
mA[0][1] = eh
for i := 1; i < n-1; i++ {
mA[i] = make([]float64, n)
mA[i][i-1] = epsilon
mA[i][i] = n2eh
mA[i][i+1] = eh
}
mA[n-1] = make([]float64, n)
mA[n-1][n-2] = epsilon
mA[n-1][n-1] = n2eh
return mA
}
func linerB(a float64) []float64 {
ah2 := a * h * h
b := make([]float64, n)
for i := 0; i < n; i++ {
b[i] = ah2
}
return b
}
func accuAns(x []float64, a float64, epsilon float64) []float64 {
xrt := make([]float64, len(x))
for i, xi := range x {
xrt[i] = xi*a + (1-a)*(1-math.Pow(math.E, (-xi/epsilon)))/
(1-math.Pow(math.E, (-1.0/epsilon)))
}
return xrt
}
func Norm1(x1, x2 []float64) (xo float64) {
for i := 0; i < len(x1); i++ {
xo += math.Abs(x1[i] - x2[i])
}
return xo
}
func Jacobi(a [][]float64, b []float64, x0 []float64, accu float64) []float64 {
size := len(b)
xrt := make([]float64, size)
x := make([]float64, size)
copy(x, x0)
var i, j int
for {
for i = 0; i < size; i++ {
xrt[i] = b[i]
for j = 0; j < i; j++ {
xrt[i] -= a[i][j] * x[j]
}
for j = i + 1; j < size; j++ {
xrt[i] -= a[i][j] * x[j]
}
xrt[i] /= a[i][i]
}
//test accuracy
if Norm1(xrt, x) < accu {
break
}
copy(x, xrt)
// fmt.Println(xrt)
}
return xrt
}
func GaussSeidel(a [][]float64, b []float64, x0 []float64, accu float64) []float64 {
size := len(b)
xrt := make([]float64, size)
x := make([]float64, size)
copy(x, x0)
var i, j int
for {
for i = 0; i < size; i++ {
xrt[i] = b[i]
for j = 0; j < i; j++ {
xrt[i] -= a[i][j] * xrt[j]
}
for j = i + 1; j < size; j++ {
xrt[i] -= a[i][j] * x[j]
}
xrt[i] /= a[i][i]
}
if Norm1(xrt, x) < accu {
break
}
copy(x, xrt)
}
return xrt
}
func SOR(a [][]float64, b []float64, omega float64, x0 []float64, accu float64) []float64 {
size := len(b)
xrt := make([]float64, size)
x := make([]float64, size)
copy(x, x0)
var i, j int
for {
for i = 0; i < size; i++ {
xrt[i] = b[i]
for j = 0; j < i; j++ {
xrt[i] -= a[i][j] * xrt[j]
}
for j = i; j < size; j++ {
xrt[i] -= a[i][j] * x[j]
}
xrt[i] = xrt[i] * omega / a[i][i]
xrt[i] += x[i]
}
if Norm1(xrt, x) < accu {
break
}
copy(x, xrt)
}
return xrt
}
func main() {
a := matrixA(1.0)
b := linerB(0.5)
x := make([]float64, n)
for i := 1; i < n; i++ {
x[i] = float64(i) / 100
}
fmt.Println(x)
fmt.Println(accuAns(x, 0.5, 1.0))
fmt.Println("---------Jacobi-----------")
fmt.Println(Jacobi(a, b, x, 1e-8))
fmt.Println("---------GaussSeidel-----------")
fmt.Println(GaussSeidel(a, b, x, 1e-8))
fmt.Println("---------SOR-----------")
fmt.Println(SOR(a, b, 1.1, x, 1e-8))
}
Raw
lab8.go
/*
* @Author: Bai Shuai
* @Date: 2014-05-20 15:52:35
* @Last Modified by: bai
* @Last Modified time: 2014-05-22 09:00:46
*/
/*
go run lab8.go 直接运行，结果如下
12.254320584751564 [0.6740198074615602 -1 0.8895596427028144]
98.52169772379699 [-0.6039723423311808 1 -0.2511351305264881 0.14895344557473505]
*/
package main
import (
"fmt"
"math"
)
var (
mA = [][]float64{
{5, -4, 1},
{-4, 6, -4},
{1, -4, 7}}
mB = [][]float64{
{25, -41, 10, -6},
{-41, 68, -17, 10},
{10, -17, 5, -3},
{-6, 10, -3, 2}}
)
// matMulLi x = Ab
func matMulLi(mta [][]float64, b, x []float64) {
length := len(b)
var tmp float64
for i := 0; i < length; i++ {
tmp = 0
for j := 0; j < length; j++ {
tmp += mta[i][j] * b[j]
}
x[i] = tmp
}
}
func macv(v []float64) float64 {
max := 0.0
for i := 0; i < len(v); i++ {
if math.Abs(v[i]) > max {
max = math.Abs(v[i])
}
}
return max
}
func lineDivF(v, u []float64, f float64) {
for i := 0; i < len(v); i++ {
u[i] = v[i] / f
}
}
// PowMethod 幂法迭代求解特征值与对应特征向量
func PowMethod(mx [][]float64, init []float64, accu float64) (lambda float64, x []float64) {
length := len(init)
k, k1 := 0, 1
v := make([][]float64, 2)
u := make([][]float64, 2)
v[0] = make([]float64, length)
v[1] = make([]float64, length)
u[0] = make([]float64, length)
u[1] = make([]float64, length)
copy(v[0], init)
copy(u[0], init)
for math.Abs(macv(v[k1])-macv(v[k])) >= accu {
//v[k1] = mx*u[k]
matMulLi(mx, u[k], v[k1])
//uk = vk/(vk)
lineDivF(v[k1], u[k1], macv(v[k1]))
k1, k = k, k1
}
return macv(v[k]), u[k]
}
func main() {
lmd, x := PowMethod(mA, []float64{1, 1, 1}, 1e-5)
fmt.Println(lmd, x)
lmd, x = PowMethod(mB, []float64{1, 1, 1, 1}, 1e-5)
fmt.Println(lmd, x)
}
Raw
lab9.go
/*
* @Author: bai
* @Date: 2014-06-02 21:07:30
* @Last Modified by: bai
* @Last Modified time: 2014-06-02 21:07:40
*/
package main
import "fmt"
func main() {
}
Raw
NA数值分析实验.md

2014年数值分析（蔡懿慈） 实验题解

运行实验代码需要安装go
安装方法

  • ubuntu sudo apt-get install golang
  • archlinux sudo pacman -S go
  • osx brew install go
  • windows 下载二进制文件安装

每个实验独立一个go文件，相互之间没有依赖
对于无输入参数的可直接运行go run labx.go执行程序,labx.go为文件名
对于需要输入运行时参数的程序，先执行go build labx.go编译代码， 再输入./labx -flag=xx，具体的参数参考对应文件的说明

截断误差与舍入误差 lab1.go

利用级数逼近在限定误差范围求ln(2)的值

多项式插值法 lab2.go 未完成

Lagrange 差值, Runge 现象
三次样条差值

曲线拟合的最小二乘法 lab3.go

使用正交多项式作最小二乘拟合，这种方法不需要求解线性方程组，只用递推公式， 并且当逼近次数增加一次时，只需要把程序中循环数加1，其余不用改变， 这是目前用多项式做曲线拟合最好的计算方法

数值积分 lab4.go

  • 复合 Simpson 求积公式
  • Romberg 外推方法求积分
  • 复合 Gauss 公式作近似积分

非线性方程求根 lab5.go

二分法、牛顿法（切线法）、迭代法
迭代函数、迭代初值对收敛性的影响

线性方程组的直接解法 lab6.go

  • 矩阵三角分解(LU 分解)的方法
  • 微小偏差对解的影响（病态矩阵）
  • 平方根法(Cholesky 分解)

线性方程组的迭代解法 lab7.go 未完成

用Jacobi法，Gauss-Seidel法和SOR法求解线性方程组的解

矩阵特征值问题 lab8.go

幂法求出矩阵按模最大的特征值与对应的特征向量

常微分方程初值问题数值解法 lab9.go 未完成

经典的四阶 Runge-Kutta 公式
四阶 Hamming 公式
Milne 公式
预测矫正系统

本实习题目会遇到两个问题:
a)按给定步长积分时,积分终点与区间端点不重合。
b)多步法改变步长。 试提出相应解决方法并实现之。

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment