📄 nihe.bas
字号:
Attribute VB_Name = "Fitting"
Option Explicit
Sub LineNiHe(x() As Double, y() As Double, a As Double, b As Double, t As Integer)
Dim i As Integer
Dim m As Double, n As Double, c As Double, d As Double
Dim sxx As Double, sxy As Double
For i = 1 To t
m = m + x(i): n = n + y(i): c = c + x(i) * x(i): d = d + y(i) * x(i)
Next i
sxy = d - n * m / t: sxx = c - m * m / t
a = sxy / sxx
b = (n - a * m) / t
End Sub
'返回值表示逼近度
'0为完全吻合,1为完全不吻合
Function LineNiHeEx(x() As Double, y() As Double, a As Double, b As Double) As Double
Dim i As Integer, t As Integer
Dim m As Double, n As Double, c As Double, d As Double, e As Double
Dim sxx As Double, sxy As Double, syy As Double
t = UBound(x)
For i = 1 To t
m = m + x(i): n = n + y(i): c = c + x(i) * x(i): d = d + y(i) * x(i): e = e + y(i) * y(i)
Next i
sxy = d - n * m / t: sxx = c - m * m / t: syy = e - n * n / t
a = sxy / sxx
b = (n - a * m) / t
LineNiHeEx = 1 - Sqr(sxy * sxy / sxx / syy)
End Function
'最小二乘法
Sub TwoMulti(a() As Double, x() As Double, n As Integer, m As Integer)
Dim i As Integer, j As Integer
Dim p() As Double, b() As Double, t As Double, s As Double, r As Double, d As Double
Dim k As Integer, z As Integer
ReDim p(m), b(n * (n + 1) / 2)
For i = 1 To n
s = 0
For j = 1 To m
t = a(j, i)
s = t * t + s
p(j) = t
Next j
d = s
For j = i + 1 To n + 1
s = 0
For z = 1 To m
s = s + p(z) * a(z, j)
Next z
r = s / d
For z = 1 To m
a(z, j) = a(z, j) - r * p(z)
Next z
k = k + 1
b(k) = r
Next j
Next i
x(n) = b(k)
For i = n - 1 To 1 Step -1
k = k - 1
s = b(k)
For j = n To i + 1 Step -1
k = k - 1
s = s - x(j) * b(k)
Next j
x(i) = s
Next i
End Sub
'建立最小二乘矩阵
't(i)=x^(i-1)为基函数
Sub CreateMinTwo(a() As Double, x() As Double, y() As Double, n As Integer, m As Integer)
Dim i As Integer, j As Integer
Dim k As Double, s As Double
For i = 1 To m
k = 1
s = x(i): a(i, n) = 1
For j = n - 1 To 1 Step -1
k = k * s
a(i, j) = k
Next j
a(i, n + 1) = y(i)
Next i
End Sub
'示例
'Private Sub tval_Click()
'Dim a(6, 4) As Double
'FillArray 6, 4, a(), 1, 0, 0, 1, 0, 1, 0, 2, 0, 0, 1, 3, -1, 1, 0, 1, 0, -1, 1, 2, -1, 0, 1, 1
'Dim b(3) As Double, x(3) As Double
'TwoMulti a(), x(), 3, 6
'Dim i As Integer
'MsgBox OutPutArray(x(), 3)
'End Sub
'三次插值第二问题(g为二阶导数)
'需给出f0-fn,c()范围-1至N+1
Sub SimpleT2(f() As Double, c() As Double, g0 As Double, gn As Double, h As Double, n As Long)
Dim i As Long
Dim y() As Double, z() As Double, a As Double
ReDim y(n), z(n)
c(0) = f(0) - h * h * g0 / 6
c(n) = f(n) - h * h * gn / 6
y(1) = 4: z(1) = f(1) * 6 - c(0)
For i = 2 To n - 1
a = 1 / y(i - 1)
y(i) = 4 - a
z(i) = f(i) * 6 - a * z(i - 1)
Next i
c(n - 1) = (z(n - 1) - c(n)) / y(n - 1)
For i = n - 2 To 1 Step -1
c(i) = (z(i) - c(i + 1)) / y(i)
Next i
c(-1) = c(0) + c(0) - c(1) + h * h * g0
c(n + 1) = c(n) + c(n) - c(n - 1) + h * h * gn
End Sub
'三次插值第一问题(g为一阶导数)
'需给出f0-fn,c()范围-1至N+1
Sub SimpleT1(f() As Double, c() As Double, g0 As Double, gn As Double, h As Double, n As Long)
Dim i As Long
Dim y() As Double, z() As Double, a As Double
ReDim y(n), z(n)
y(0) = 4: z(0) = f(0) * 6 + (h + h) * g0: y(1) = 3.5: z(1) = f(1) * 6 - z(0) / 4
For i = 2 To n - 1
a = 1 / y(i - 1)
y(i) = 4 - a
z(i) = f(i) * 6 - a * z(i - 1)
Next i
a = 2 / y(n - 1)
c(n) = (f(n) * 6 - (h + h) * gn - a * z(n - 1)) / (4 - a)
For i = n - 1 To 1 Step -1
c(i) = (z(i) - c(i + 1)) / y(i)
Next i
c(0) = (z(0) - c(1) - c(1)) / y(0)
c(-1) = c(1) - (h + h) * g0
c(n + 1) = c(n - 1) + (h + h) * gn
End Sub
'三次插值第三问题
'周期函数
'需给出f0-fn,c()范围-1至N+1
Sub SimpleT3(f() As Double, c() As Double, h As Double, n As Long)
Dim m() As Double, k() As Double, s() As Double, t() As Double
ReDim m(n), k(n), s(n), t(n)
Dim i As Long, l As Double, s1 As Double, s2 As Double
m(1) = 4: k(1) = 1 / 4: s(1) = 1: t(1) = 6 * f(1)
For i = 2 To n - 1
l = 1 / m(i - 1)
m(i) = 4 - l: k(i) = -k(i - 1) / m(i)
s(i) = -l * s(i - 1)
t(i) = f(i) * 6 - l * t(i - 1)
Next i
l = 1 / m(n - 1) + k(n - 1)
For i = 1 To n - 2
s1 = k(i) * s(i) + s1: s2 = k(i) * t(i) + s2
Next i
c(n) = (f(n) * 6 - s2 - l * t(n - 1)) / (4 - s1 - l * (1 + s(n - 1)))
l = c(n): c(n - 1) = (t(n - 1) - (1 + s(n - 1)) * l) / m(n - 1)
For i = n - 2 To 1 Step -1
c(i) = (t(i) - l * s(i) - c(i + 1)) / m(i)
Next i
c(0) = c(n): c(-1) = c(n - 1): c(n + 1) = c(1)
End Sub
'调用示例
'Private Sub tval_Click()
'Dim i As Integer
'Dim f(18) As Double, c(-1 To 19) As Double, t As Double
'For i = 0 To 18
'f(i) = Sin(PI / 18 * i)
'Next i
'SimpleT3 f(), c(), PI / 18, 18
'For i = -1 To 19
't = t + c(i) * omi3(36 / PI - i)
'Next i
'tval = t
'End Sub
'获取切比雪夫多项式系数a(k)
Function QieBi(k As Integer, n As Long) As Double
Dim i As Long
Dim t As Double, x As Double, s As Double
For i = 1 To n
t = (i + i - 1) * PI / (n + n)
x = Cos(t)
s = s + fqie(x) * Cos(k * t)
Next i
QieBi = (s + s) / n
If k = 0 Then QieBi = QieBi / 2
End Function
Function fqie(x As Double) As Double
fqie = Cos(x * PI / 4)
End Function
Function Pade(c() As Double, a() As Double, b() As Double, m As Integer, n As Integer) As Boolean
Dim i As Integer, j As Integer, q As Integer
Dim y() As Double, d() As Double, t As Double
ReDim d(n, n), y(n)
For i = 1 To n
If m + i > n Then q = 1 Else q = n + 1 - m - i
For j = q To n
d(i, j) = c(i + j + m - n - 1)
Next j
y(i) = -c(m + i)
Next i
If EquGaussEx(d(), y(), b(), n) = False Then Exit Function
For i = 1 To n \ 2
t = b(i)
b(i) = b(n + 1 - i)
b(n + 1 - i) = t
Next i
For i = 0 To m
If i > n Then q = n Else q = i
t = 0
For j = 1 To q
t = t + b(j) * c(i - j)
Next j
a(i) = t + c(i)
Next i
b(0) = 1
Pade = True
End Function
'例程
'Private Sub tval_Click()
'Dim i As Integer
'Dim c(5) As Double, a(3) As Double, b(2) As Double
'c(0) = 0: c(1) = 1: c(2) = -0.5: c(3) = 1 / 3: c(4) = -0.25: c(5) = 0.2
'Pade c(), a(), b(), 3, 2
'For i = 0 To 3
'tval = tval & a(i) & vbcrlf
'Next i
'For i = 0 To 2
'tval = tval & b(i) & vbcrlf
'Next i
'End Sub
'需提供个m+2n参数
'尚未调试
Function PadeEx(c() As Double, a() As Double, b() As Double, m As Integer, n As Integer) As Boolean
Dim i As Integer, j As Integer, k As Integer
Dim y() As Double, d() As Double, t As Double
For i = 1 To n
k = m + i
For j = 1 To n
d(i, j) = c(j + k) + c(Abs(k - j))
Next j
If k <= n Then d(i, k) = d(i, k) + c(0)
y(i) = -c(k) - c(k)
Next i
If EquGaussEx(d(), y(), b(), n) = False Then Exit Function
b(0) = 1
For i = 1 To m
t = 0
For j = 1 To n
t = t + (c(i + j) + c(Abs(i - j))) * b(j)
Next j
If i <= n Then t = t + b(i) * c(0)
a(i) = c(i) + t / 2
Next i
t = 0
For i = 1 To n
t = t + b(i) + c(i)
Next i
a(0) = c(0) + t / 2
PadeEx = True
End Function
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -