79 lines
2.3 KiB
Plaintext
79 lines
2.3 KiB
Plaintext
Sub GaussJordan(matrix() As Double,rhs() As Double,ans() As Double)
|
|
Dim As Integer n=Ubound(matrix,1)
|
|
Redim ans(0):Redim ans(1 To n)
|
|
Dim As Double b(1 To n,1 To n),r(1 To n)
|
|
For c As Integer=1 To n 'take copies
|
|
r(c)=rhs(c)
|
|
For d As Integer=1 To n
|
|
b(c,d)=matrix(c,d)
|
|
Next d
|
|
Next c
|
|
#macro pivot(num)
|
|
For p1 As Integer = num To n - 1
|
|
For p2 As Integer = p1 + 1 To n
|
|
If Abs(b(p1,num))<Abs(b(p2,num)) Then
|
|
Swap r(p1),r(p2)
|
|
For g As Integer=1 To n
|
|
Swap b(p1,g),b(p2,g)
|
|
Next g
|
|
End If
|
|
Next p2
|
|
Next p1
|
|
#endmacro
|
|
For k As Integer=1 To n-1
|
|
pivot(k) 'full pivoting
|
|
For row As Integer =k To n-1
|
|
If b(row+1,k)=0 Then Exit For
|
|
Var f=b(k,k)/b(row+1,k)
|
|
r(row+1)=r(row+1)*f-r(k)
|
|
For g As Integer=1 To n
|
|
b((row+1),g)=b((row+1),g)*f-b(k,g)
|
|
Next g
|
|
Next row
|
|
Next k
|
|
'back substitute
|
|
For z As Integer=n To 1 Step -1
|
|
ans(z)=r(z)/b(z,z)
|
|
For j As Integer = n To z+1 Step -1
|
|
ans(z)=ans(z)-(b(z,j)*ans(j)/b(z,z))
|
|
Next j
|
|
Next z
|
|
End Sub
|
|
|
|
'Interpolate through points.
|
|
Sub Interpolate(x_values() As Double,y_values() As Double,p() As Double)
|
|
Var n=Ubound(x_values)
|
|
Redim p(0):Redim p(1 To n)
|
|
Dim As Double matrix(1 To n,1 To n),rhs(1 To n)
|
|
For a As Integer=1 To n
|
|
rhs(a)=y_values(a)
|
|
For b As Integer=1 To n
|
|
matrix(a,b)=x_values(a)^(b-1)
|
|
Next b
|
|
Next a
|
|
'Solve the linear equations
|
|
GaussJordan(matrix(),rhs(),p())
|
|
End Sub
|
|
|
|
'======================== SET UP THE POINTS ===============
|
|
|
|
Dim As Double x(1 To ...)={0,1,2,3,4,5,6,7,8,9,10}
|
|
Dim As Double y(1 To ...)={1,6,17,34,57,86,121,162,209,262,321}
|
|
|
|
Redim As Double Poly(0)
|
|
'Get the polynomial Poly()
|
|
Interpolate(x(),y(),Poly())
|
|
|
|
'print coefficients to console
|
|
print "Polynomial Coefficients:"
|
|
print
|
|
For z As Integer=1 To Ubound(Poly)
|
|
If z=1 Then
|
|
Print "constant term ";tab(20);Poly(z)
|
|
Else
|
|
Print tab(8); "x^";z-1;" = ";tab(20);Poly(z)
|
|
End If
|
|
Next z
|
|
|
|
sleep
|