d2r3.frm

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

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   5400
   ClientLeft      =   2055
   ClientTop       =   750
   ClientWidth     =   4725
   LinkTopic       =   "Form1"
   ScaleHeight     =   5400
   ScaleWidth      =   4725
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3240
      TabIndex        =   0
      Top             =   4800
      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 D2R3
    'DRIVER for routine SPLINE
    N = 20
    PI = 3.1415926
    Dim X(20), Y(20), Y2(20)
    Print
    Print Tab(5); "Second-derivative for sin(x) from 0 to PI"
    '生成插值数组
    For I = 1 To 20
        X(I) = I * PI / N
        Y(I) = Sin(X(I))
    Next I
    '用三次样条函数计算二阶导数
    YP1 = Cos(X(1))
    YPN = Cos(X(N))
    Call SPLINE(X(), Y(), N, YP1, YPN, Y2())
    '验证计算结果
    Print Tab(5); "            spline        actual"
    Print Tab(5); " angle     2nd deriv     2nd deriv"
    For I = 1 To N
        Print Tab(6); Format$(X(I), "##.##");
        Print Tab(16); Format$(Y2(I), "##.#####0");
        Print Tab(30); Format$(-Sin(X(I)), "##.#####0")
    Next I
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