d2r4.frm

VB常用数值算法集 内含有解线性代数方程组

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   6240
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   4680
   LinkTopic       =   "Form1"
   ScaleHeight     =   6240
   ScaleWidth      =   4680
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3120
      TabIndex        =   0
      Top             =   5520
      Width           =   1215
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Private Sub Command1_Click()
    'PROGRAM D2R4
    'DRIVER for routine SPLINT,which calls SPLINE
    NP = 10
    PI = 3.141593
    Dim XA(10), YA(10), Y2(10)
    For NFUNC = 1 To 2
        If NFUNC = 1 Then
            Print
            Print Tab(5); "Sine function from 0 to PI"
            For I = 1 To NP
                XA(I) = I * PI / NP
                YA(I) = Sin(XA(I))
            Next I
            YP1 = Cos(XA(1))
            YPN = Cos(XA(NP))
        ElseIf NFUNC = 2 Then
            Print
            Print Tab(5); "Exponential function from 0 to 1"
            For I = 1 To NP
                XA(I) = 1# * I / NP
                YA(I) = Exp(XA(I))
            Next I
            YP1 = Exp(XA(1))
            YPN = Exp(XA(NP))
        Else
            Stop
        End If
        'Call SPLINE to get second derivatives
        Call SPLINE(XA(), YA(), NP, YP1, YPN, Y2())
        'Call SPLINT for interpolations
        Print Tab(5); "   x            f(x)     interpolation"
        For I = 1 To 10
            If NFUNC = 1 Then
                X = (-0.05 + I / 10#) * PI
                F = Sin(X)
            ElseIf NFUNC = 2 Then
                X = (-0.05 + I / 10#)
                F = Exp(X)
            End If
            Call SPLINT(XA(), YA(), Y2(), NP, X, Y)
            Print Tab(7); Format$(X, "##.##0000");
            Print Tab(19); Format$(F, "##.###000");
            Print Tab(32); Format$(Y, "##.###000")
        Next I
        Print Tab(5); "**********************************"
    Next NFUNC
End Sub
   Sub SPLINE(X(), Y(), N, YP1, YPN, Y2())
      Dim U(100)
      If YP1 > 9.9E+29 Then
          Y2(1) = 0
          U(1) = 0
      Else
          Y2(1) = -0.5
          AAA = (Y(2) - Y(1)) / (X(2) - X(1))
          U(1) = (3# / (X(2) - X(1))) * (AAA - YP1)
      End If
      For I = 2 To N - 1
          SIG = (X(I) - X(I - 1)) / (X(I + 1) - X(I - 1))
          P = SIG * Y2(I - 1) + 2#
          Y2(I) = (SIG - 1#) / P
          AAA = (Y(I + 1) - Y(I)) / (X(I + 1) - X(I))
          BBB = (Y(I) - Y(I - 1)) / (X(I) - X(I - 1))
          CCC = X(I + 1) - X(I - 1)
          U(I) = (6# * (AAA - BBB) / CCC - SIG * U(I - 1)) / P
      Next I
      If YPN > 9.9E+29 Then
          QN = 0#
          UN = 0#
      Else
          QN = 0.5
          AAA = YPN - (Y(N) - Y(N - 1)) / (X(N) - X(N - 1))
          UN = (3# / (X(N) - X(N - 1))) * AAA
      End If
      Y2(N) = (UN - QN * U(N - 1)) / (QN * Y2(N - 1) + 1#)
      For K = N - 1 To 1 Step -1
          Y2(K) = Y2(K) * Y2(K + 1) + U(K)
      Next K
   End Sub
Sub SPLINT(XA(), YA(), Y2A(), N, X, Y)
    KLO = 1
    KHI = N
1   If KHI - KLO > 1 Then
        K = (KHI + KLO) / 2
        If XA(K) > X Then
            KHI = K
        Else
            KLO = K
        End If
        GoTo 1
    End If
    H = XA(KHI) - XA(KLO)
    If H = 0 Then
        Print "  PAUSE  'BAD  XA  INPUT'"
        Exit Sub
    End If
    A = (XA(KHI) - X) / H
    B = (X - XA(KLO)) / H
    AAA = A * YA(KLO) + B * YA(KHI)
    BBB = (A ^ 3 - A) * Y2A(KLO) + (B ^ 3 - B) * Y2A(KHI)
    Y = AAA + BBB * (H ^ 2) / 6#
End Sub