有偿求“最小二乘法曲线拟合”计算软件

没有这方面的经历，在网上找了段代码，可以运行，拟合的是二次曲线。
新建一工程，窗口上放4个按钮，从command1到command4分别是“画坐标、采样点生成、曲线拟合、退出”，复制下列代码并运行即可看到效果。
Option Explicit
Dim x()  As Double, y() As Double
Dim A(20, 20) As Double, M As Double, B() As Double  '最多取20次的拟合
Dim N As Double, I As Double, j As Double
Dim xiaoA() As Double
Dim Xmin As Double, Xmax As Double
Dim Ymin As Double, Ymax As Double
Dim X0pos As Double, Y0pos As Double
Dim xmaxpos As Double, ymaxpos As Double
Dim xstep As Double, ystep As Double
Dim xl As Double, yl As Double
Dim xbc As Double, ybc As Double
Dim bc As Double
Dim Xh As Double

Private Sub HuaZuoBiao(x() As Double, y() As Double)

ReDim xpos(I) As Double
ReDim ypos(I) As Double
ReDim x(I), y(I)

X0pos = Width * 0.25 '坐标原点最左点
Y0pos = Height * 0.75 '坐标原点最低点

xmaxpos = Width * 0.85 '坐标最右点
ymaxpos = Height * 0.15 '坐标最高点

xstep = (xmaxpos - X0pos) / (Xmax - Xmin) '对应X轴上单位长度代表的屏幕宽度值
ystep = (ymaxpos - Y0pos) / (Ymax - Ymin) '对应Y轴上单位长度代表的屏幕高度值

'在屏幕上画直角坐标系

ForeColor = vbBlue
Line (Width * 0.1, Y0pos)-(Width * 0.9, Y0pos) '画X坐标轴,从左10%，到右的90%处
Line (X0pos, Height * 0.1)-(X0pos, Height * 0.9) '画y坐标轴,从上10%，到下的90%处

Font.Size = 20 '指定X轴，Y轴标志的字体大小

CurrentX = Width * 0.9
CurrentY = Y0pos + 100
Print "X" '在横线上画X轴标志
 
 '在横线上画X轴箭头标志
CurrentX = Width * 0.9
CurrentY = Y0pos
Line (CurrentX - 200, CurrentY - 50)-(CurrentX, CurrentY)
Line (CurrentX, CurrentY)-(CurrentX - 200, CurrentY + 50)

CurrentX = X0pos - 500
CurrentY = Height * 0.1
Print "y" '在纵线上画Y轴标志

 '在纵线上画Y轴箭头标志
CurrentX = X0pos
CurrentY = Height * 0.1
Line (CurrentX - 50, CurrentY + 200)-(CurrentX, CurrentY)
Line (CurrentX, CurrentY)-(CurrentX + 50, CurrentY + 200)

CurrentX = X0pos + 200 '此为Y轴左边500绝对坐标处
CurrentY = Y0pos + 400 '取当前Y轴上的相对坐标值
Print "f=f(x)" '在Y轴左边500绝对坐标处对应显示Y轴相对坐标刻度值

xl = Xmax - Xmin
yl = Ymax - Ymin

If xl < 0.01 Then
 xbc = 0.001
ElseIf xl <= 0.1 Then
 xbc = 0.01
ElseIf xl <= 2 Then
 xbc = 0.1
ElseIf xl <= 20 Then
 xbc = 1
ElseIf xl <= 120 Then
 xbc = 10
ElseIf xl <= 1000 Then
 xbc = 100
ElseIf xl <= 10000 Then
 xbc = 1000
Else
 xbc = 10000
End If

If yl < 0.01 Then
 ybc = 0.001
ElseIf yl <= 0.1 Then
 ybc = 0.01
ElseIf yl <= 2 Then
 ybc = 0.1
ElseIf yl <= 20 Then
 ybc = 1
ElseIf yl <= 120 Then
 ybc = 10
ElseIf yl <= 1000 Then
 ybc = 100
ElseIf yl <= 10000 Then
 ybc = 1000
Else
 ybc = 10000
End If
  

For bc = Xmin To Xmax Step xbc

If bc <= Xmax Then
x(j) = bc 'X轴上的相对坐标值

xpos(j) = X0pos + (x(j) - Xmin) * xstep

Line (xpos(j), Y0pos)-(xpos(j), ymaxpos), vbRed ' 画垂直于X轴的刻度线，只画了100个绝对尺寸

Else
End If

Font.Size = 10 '指定X轴，Y轴坐标刻度值的字体大小

CurrentX = xpos(j) - 200 '取当前X轴上的相对坐标值
CurrentY = Y0pos + 100 '此为X轴下方100绝对坐标处
Print x(j) '在X轴下方100绝对坐标处对应显示X轴相对坐标刻度值

Next bc

For bc = Ymin To Ymax Step ybc

If bc <= Ymax Then
y(j) = bc 'X轴上的相对坐标值

ypos(j) = Y0pos + (y(j) - Ymin) * ystep

Line (X0pos, ypos(j))-(xmaxpos, ypos(j)), vbRed ' 画垂直于X轴的刻度线，只画了100个绝对尺寸

Else
End If

Font.Size = 10 '指定X轴，Y轴坐标刻度值的字体大小

CurrentX = X0pos - 500 '取当前X轴上的相对坐标值
CurrentY = ypos(j) - 100 '此为X轴下方100绝对坐标处
Print y(j) '在X轴下方100绝对坐标处对应显示X轴相对坐标刻度值

Next bc

End Sub

Private Sub ZuoDian(x() As Double, y() As Double)


ReDim xpos(I) As Double
ReDim ypos(I) As Double

For I = 0 To N


xpos(I) = X0pos + (x(I) - Xmin) * xstep
ypos(I) = Y0pos + (y(I) - Ymin) * ystep

If y(I) <= Ymax Then

   DrawWidth = 4
   PSet (xpos(I), ypos(I)), vbRed

Else
End If
  
Next I

  DrawWidth = 1

End Sub

Private Sub HuaQuXian(xiaoA() As Double)

ReDim xpos(I) As Double
ReDim ypos(I) As Double

Dim Ysum As Double, Ii As Double
For Ii = Xmin To Xmax Step 1 / (Xmax - Xmin)

      Ysum = 0
      For j = 1 To M
          Ysum = Ysum + xiaoA(j) * Ii ^ (j - 1)
      Next j
      
          xpos(I) = X0pos + (Ii - Xmin) * xstep
          ypos(I) = Y0pos + (Ysum - Ymin) * ystep
DrawWidth = 2
If Ii = Xmin Then
          xpos(0) = X0pos + (Ii - Xmin) * xstep
          ypos(0) = Y0pos + (Ysum - Ymin) * ystep
PSet (xpos(0), ypos(0))
Else
End If

If Ysum <= Ymax Then

DrawWidth = 2
Line -(xpos(I), ypos(I)), vbBlue

Else
End If

Next Ii
DrawWidth = 1
     
End Sub
Private Sub JieFangCheng(A() As Double, B() As Double, x() As Double)
Dim nn As Double
nn = UBound(B)

Dim TempA As Double, L As Double, K As Double, Kk As Double
Dim Ii As Double, ChuShu As Double, Sum As Double
For I = 1 To nn
    L = 0: Kk = 0
    For j = I To nn
      If A(j, I) = 0 Then L = L + 1
    Next j
    For j = I To nn - L
      If A(j, I) = 0 Then
        Kk = Kk + 1
        For K = I To nn
           TempA = A(j, K)
           A(j, K) = A(nn - Kk + 1, K)
           A(nn - Kk + 1, K) = TempA
        Next K
        TempA = B(j): B(j) = B(nn - Kk + 1): B(nn - Kk + 1) = TempA
      End If
    Next j
              
    For Ii = I To nn - L
      ChuShu = A(Ii, I)
      For j = I To nn
         A(Ii, j) = A(Ii, j) / ChuShu
      Next j
      B(Ii) = B(Ii) / ChuShu
    Next Ii
    For Ii = I + 1 To nn - L
      For j = I To nn
         A(Ii, j) = A(Ii, j) - A(I, j)
      Next j
      B(Ii) = B(Ii) - B(I)
    Next Ii
Next I
For I = 1 To nn
    For j = 1 To I - 1
      A(I, j) = 0
    Next j
Next I
      
      
x(nn) = B(nn) / A(nn, nn)
For I = nn - 1 To 1 Step -1
   Sum = 0
   For j = I + 1 To nn
      Sum = Sum + A(I, j) * x(j)
   Next j
   x(I) = (B(I) - Sum) / A(I, I)
Next I

End Sub

Private Sub Command1_Click()
Cls

Xmin = 0  ' InputBox("请输入x坐标下限值", "x坐标下限值", 0)
Ymin = 0  'InputBox("请输入y坐标下限值", "y坐标下限值", 0)
Xmax = 10 ' InputBox("请输入x坐标上限值", "x坐标上限值度", 10)
Ymax = 10 'InputBox("请输入y坐标上限值", "y坐标上限值度", 10)
N = 20
For I = 0 To N

 ReDim Preserve x(I)
 ReDim Preserve y(I)
     
Next I
Call HuaZuoBiao(x, y)

End Sub

Private Sub Command2_Click()

For I = 0 To N
  
 x(I) = Xmin + I * (Xmax - Xmin) / N 'InputBox("请输入X坐标测量值", "X坐标值", "0") '
     
 y(I) = Sin(x(I)) + 5 ' InputBox("请输入Y坐标测量值", "Y坐标值", "0") '
      
Next I

Call ZuoDian(x, y)


End Sub

Private Sub Command3_Click()

M = 20 'InputBox("请输入拟合曲线次数M", "拟合曲线", 3)
Erase B: Erase xiaoA: Erase A   '必不可少***********
ReDim B(M): ReDim xiaoA(1 To M)
'形成方程组的各元素
A(1, 1) = N
For I = 1 To N
   B(1) = B(1) + y(I)
Next I
For j = 2 To M
   For I = 1 To N
      A(1, j) = A(1, j) + x(I) ^ (j - 1)
   Next I
Next j
For I = 2 To M
   For j = 1 To M
      For Xh = 1 To N
         A(I, j) = A(I, j) + x(Xh) ^ (I + j - 2)
         If j = 1 Then
            B(I) = B(I) + x(Xh) ^ (I - 1) * y(Xh)
         End If
      Next Xh
   Next j
Next I

Call JieFangCheng(A, B, xiaoA)

ForeColor = vbBlack
PSet (0, 0)
For I = 1 To M
   'Print Tab(6); "a"; I - 1; Tab(12); "="; xiaoA(I);
Next I
Dim Str As String: Str = "y="
For I = 1 To M    '写方程
   If I < M Then
       Str = Str & xiaoA(I) & "*x^" & I - 1 & "+"
   Else
       Str = Str & xiaoA(I) & "*x^" & I - 1
   End If
Next I
Print vbCrLf; "曲线方程:"; vbCrLf & Str

Call HuaQuXian(xiaoA)

End Sub

Private Sub Command4_Click()
End
End Sub

Private Sub Form_Load()

Width = Screen.Width * 1 '取屏幕宽度的一半
'Height = Screen.Height * 0.5 '取屏幕高度的一半
Height = Screen.Width * 1 '取屏幕宽度的一半
Left = (Screen.Width - Width) / 2 '使窗体居屏幕中心
Top = (Screen.Height - Height) / 2 '使窗体居屏幕中心

End Sub