|
|
|
|
2018年12月02日 |
|
隱函數求解
如果您需要計算一個類似x^3+2*x*y^2+x*y^3-x^2*z^2-25=0之多變數函數在變數{x,-5~5},{y,-5~5},{z,-5~5}範圍時之各種解或曲線，如果您手上沒有任何電腦軟體時？該怎麼辦？在此我介紹一個非常簡便的方法？如果您會多項式單變數解法，那麼問題就解決一半。以一個三變數x，y，z之多項式來說，我們可以先假定z值與y值後，解單變數多項式即可。至於單變數多項式電腦解法在網路上可以搜索到(C++、VB、或JAVA都有)。在本報告中，我們可以利用到Vbscript控制項，讓User直接輸入公式，如 ” x^3+2*x*y^2+x*y^3-x^2*z^2-25 ”，然後其他計算工作：函數次數(order)，多項式係數，求根等，就交由電話代勞。輸入時，公式應符合VB程式語法外，x變數之冪次規定比較嚴格，如x三次方，應寫成x^3而非x*x*x，但y、z則可以任意，故y^2與y*y同義，x*z^2與x*z*z或x^1*z^2相當。下面就是部份電腦程式。
Private Sub FindPolyCoeffLH(eqStIn As String, yIn As Single, zIn As Single, eqStOut As String, nPoly As Integer, Acoeff_0() As Double, Optional iTest As Boolean = False) '*********************************************** 'Please note here function :A(0)x^n+A(1)x^(n-1)+A(2)x^(n-2)+........+A(n) 'if 2x^5+3x^4+2x^3+3x^2+4x+5=0 : then A(0)=2,A(1)=3,A(2)=2,A(3)=3,A(4)=4,A(5)=5 '*************************************************** Dim i As Integer, count As Integer, j As Integer Dim x As Single, y As Single, z As Integer Dim ret As Variant '----------------------------------------------- Dim st As String st = eqStIn ret = FunValue(st, 0, 0, 0) 'List1.AddItem "a(0)=" & ret st = LCase(st) st = Replace(st, "x*", "x^1*") stArray = Split(st, "+", -1, vbTextCompare) 'List1.AddItem St count = 0 For i = 0 To UBound(stArray) ReDim Preserve Axyz(0 To count) Axyz(count) = stArray(i) 'List1.AddItem "i= " & count & " ;Axyz(" & count & ")= " & Axyz(count) count = count + 1 Next i '------------------------ count = count - 1 nPoly = -999
For i = 0 To count For j = 1 To 10 If InStr(Axyz(i), Trim("x^" & j)) And j >= nPoly Then nPoly = j End If Next j Next i 'List1.AddItem " order of poly= " & nPoly
ReDim Preserve Acoeff_0(0 To nPoly) Acoeff_0(nPoly) = FunValue(eqStIn, 0, yIn, zIn) For j = 1 To nPoly Acoeff_0(nPoly - j) = 0 For i = 0 To count ' MsgBox ("x^" & Trim("x^" & j) & ";axyz= " & Axyz(i)) If InStr(Axyz(i), Trim("x^" & j)) Then Acoeff_0(nPoly - j) = Acoeff_0(nPoly - j) + FunValue(Axyz(i), 1, yIn, zIn) End If Next i Next j
If iTest Then For i = 0 To nPoly List1.AddItem "z= " & zIn & " ;y= " & yIn & " A(" & i & ")= " & Acoeff_0(i) Next End If Dim Fxyz As String Fxyz = Str(Acoeff_0(0)) For i = 1 To nPoly Fxyz = Fxyz & Trim("+" & Str(Acoeff_0(i)) & Trim("*x^" & i)) Next i If iTest Then List1.AddItem Fxyz eqStOut = Fxyz End Sub
The figure shown as below is the result using the snippet code attached hereafter.
Private Sub ImplicitB(EqIn As String) '*********************************************** 'Please note here function :A(0)x^n+A(1)x^(n-1)+A(2)x^(n-2)+........+A(n) 'if 2x^5+3x^4+2x^3+3x^2+4x+5=0 : then A(0)=2,A(1)=3,A(2)=2,A(3)=3,A(4)=4,A(5)=5 '*************************************************** EyeR = 10 EyeTheta = pi * 0.2 EyePhi = pi * 0.2 m3PProject glPST, m3Perspective, EyeR, EyePhi, EyeTheta, FocusX, FocusY, FocusZ, 0, 1, 0 FactorMaker glPST Dim nPoly As Integer, Acoeff_0() As Double, i As Long, nans As Integer, Nroots As Integer Dim EqOut As String, Roots As PolyRoot, xans As Single, count As Long Dim xp As Single, yp As Single, ptsIn() As pointApi, ptsOut() As pointApi Dim ptTpt As cadPoint Dim rgbR As Byte, rgbG As Byte, rgbB As Byte, yFctUse As Single Dim y As Single, z As Single, countTL As Integer Dim stColl As New Collection, stArray() As String Dim xFun() As Single, yFun() As Single, zFun() As Single Dim ansReal() As Double, ansImag() As Double Dim time1 As Double time1 = Timer EqIn = LCase(EqIn) If InStr(EqIn, "^") = 0 Then EqIn = Replace(EqIn, "x", "x^1") End If List1.AddItem EqIn countTL = 0 yFctUse = -1 Set stColl = Nothing For z = 5 To -5 Step -0.25 count = 0 Erase ptsIn, ptsOut List1.AddItem "z= " & z For y = 5 To -5 Step -0.25 Call FindPolyCoeffLH(EqIn, y, z, EqOut, nPoly, Acoeff_0) ReDim Preserve Acoeff_0(0 To nPoly) ReDim Preserve ansReal(1 To nPoly) ReDim Preserve ansImag(1 To nPoly) Call SolpolyN(nPoly, Acoeff_0, 0.0001, 0, 0, ansReal, ansImag) For i = 1 To nPoly If Abs(Round(ansImag(i), 5)) <= 0.0001 Then xans = Round(ansReal(i), 5) List1.AddItem "x,y= " & xans & " ; " & y xp = xans * glAx + y * glAy + z * glAz yp = xans * glBx + y * glBy + z * glBz Picture1.FillColor = QBColor(countTL Mod 16) Picture1.FillStyle = vbSolid Picture1.Circle (xp, yp), (0.1 + (z) * 0.0125), vbRed 'QBColor(count Mod 16) ReDim Preserve ptsIn(0 To count) ptsIn(count).x = Picture2.ScaleWidth \ 2 + CLng(Picture2.ScaleX(xp * 25 - Picture2.ScaleLeft, Picture2.ScaleMode, vbPixels)) ptsIn(count).y = Picture2.ScaleHeight \ 2 + CLng(Picture2.ScaleY(yp * 25 * yFctUse - Picture2.ScaleTop, Picture2.ScaleMode, vbPixels)) stColl.Add Trim(Round(xans, 3) & "," & Round(y, 3) & "," & Round(z, 3)) count = count + 1 End If Next i Next y
count = count - 1 If count >= 3 Then Dim nansHull As Long nansHull = MakeConvexHull(Picture2, ptsIn, ptsOut, QBColor(count Mod 16)) End If countTL = countTL + 1 List2.AddItem "*******z= " & z & " ;npt= " & count Next z SortNum1_Duplicate stColl count = 0 For i = 1 To stColl.count stArray = Split(stColl(i), ",", -1, vbTextCompare) ReDim Preserve xFun(0 To count), yFun(0 To count), zFun(0 To count) xFun(count) = Val(stArray(0)) yFun(count) = Val(stArray(1)) zFun(count) = Val(stArray(2)) count = count + 1 Next i count = count - 1 Dim FileNo As Integer, Fname As String Fname = "f:\implicitfunction\ImplicitData.dat" FileNo = FreeFile Open Fname For Output As #FileNo For i = 1 To count Write #FileNo, xFun(i), yFun(i), zFun(i) Next i Close #FileNo Lbltime.Caption = (Timer - time1) & "sec" End Sub
Sub SolpolyN(Nc As Integer, A() As Double, Tol As Double, Rzero As Double, Szero As Double, ansReal() As Double, ansImag() As Double) Dim i As Integer Dim m As Integer, j As Integer, L As Integer, n As Integer, r As Double, S As Double Dim Radtrm As Double, b(9) As Double, Rc As Double, Sc As Double Dim p(9) As Double, K As Integer, Rcr As Double, Scr As Double, Q(9) As Double, Rcs As Double Dim Scs As Double, Denom As Double, Rnum As Double, Snum As Double, Delr As Double, Dels As Double Dim NewN As Integer, rad As Double, Fa As Double, Fb As Double, Ftpt As Double 'Please note here function :A(0)x^n+A(1)X^(n-1)+A(2)x^(n-2)+........+A(n) 'if 2x^5+3x^4+2x^3+3x^2+4x+5=0 : then A(0)=2,A(1)=3,A(2)=2,A(3)=3,A(4)=4,A(5)=5 'MsgBox ("nc=" & Nc) On Error Resume Next For i = 1 To Nc A(i) = A(i) / A(0) 'MsgBox ("A(I)= " & i & ";" & A(i)) Next i A(0) = 1 m = 0 j = 1 '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Do r = Rzero 'Rzero initial S = Szero 'Szero initial L = 0 n = Nc - 2 * m ' MsgBox ("N=" & n) Select Case n '///////////////////////////////// Case Is < 2 ansReal(j) = -A(1) ' Exit Sub 'GoTo 500 Case Is = 2 Fa = A(1) * A(1) - 4# * A(2) Select Case Fa Case Is < 0# Radtrm = -(A(1) * A(1) - 4# * A(2)) rad = Sqr(Radtrm) ansReal(j) = A(1) / 2# ansReal(j + 1) = -A(1) / 2# ansImag(j) = rad / 2# ansImag(j + 1) = -rad / 2# Exit Sub 'GoTo 500 Case Is = 0# ansReal(j) = A(1) / 2# ansReal(j + 1) = -A(1) / 2# ansImag(j) = 0# ansImag(j + 1) = 0# Case Is > 0# rad = Sqr(A(1) * A(1) - 4# * A(2)) ansReal(j) = (-A(1) + rad) / 2# ansReal(j + 1) = (-A(1) - rad) / 2# ansImag(j) = 0# ansImag(j + 1) = 0# Exit Sub End Select 'calculate approx. values of R & S計算R,S近似值 Case Is > 2 End Select '//////////////////////////////////////////////////// Do b(1) = A(1) - r b(2) = A(2) - r * b(1) - S For K = 3 To n b(K) = A(K) - r * b(K - 1) - S * b(K - 2) Next K Rc = b(n - 1) Sc = b(n) + r * b(n - 1) p(1) = -1# p(2) = r - b(1) For K = 3 To n p(K) = -b(K - 1) - r * p(K - 1) - S * p(K - 2) Next K Rcr = p(n - 1) Scr = p(n) + r * p(n - 1) + b(n - 1) 'Calculate b(k) 計算 b(k)偏微分 Q(1) = 0# Q(2) = -1# For K = 3 To n Q(K) = -b(K - 2) - r * Q(K - 1) - S * Q(K - 2) Next K Rcs = Q(n - 1) Scs = Q(n) + r * Q(n - 1) ' modify Denom, Rnum, Snum計算微分修正程式 Denom = Rcr * Scs - Rcs * Scr Rnum = -Rc * Scs + Sc * Rcs Snum = -Rcr * Sc + Scr * Rc Delr = Rnum / Denom Dels = Snum / Denom '計算R,S下一階近似值 r = r + Delr S = S + Dels '測試微分修正程式是否收斂 If (Abs(Delr) - Tol) <= 0# And (Abs(Dels) - Tol) <= 0# Then Exit Do '測試迭代次數是否超過預設值 If L > 100 Then Exit Sub L = L + 1 Loop While True
'計算二次方程式F(M)之一對根值大小 Ftpt = (r * r - 4# * S) Select Case Ftpt Case Is < 0# Radtrm = -Ftpt rad = Sqr(Radtrm) ansReal(j) = -r / 2# ansReal(j + 1) = -r / 2# ansImag(j) = rad / 2# ansImag(j + 1) = -rad / 2#
Case Is = 0# ansReal(j) = -r / 2# ansReal(j + 1) = -r / 2# ansImag(j) = 0# ansImag(j + 1) = 0# Case Is > 0# rad = Sqr(Ftpt) ansReal(j) = (-r + rad) / 2# ansReal(j + 1) = (-r - rad) / 2# ansImag(j) = 0# ansImag(j + 1) = 0# End Select 'P(N-2)取代P(N)後繼續求解 m = m + 1 j = j + 2 NewN = Nc - 2 * m For K = 1 To NewN A(K) = b(K) Next K Loop While True '!!!!!!!!!!!!!!!!!!!!!!!! End Sub
下面是利用筆者所寫的程式所計算之成果畫面，左邊為利用點(浮點資料)資料所畫出者，右邊則進一步以點資料，利用凸多邊形MakeConvexHull所畫出者。
此因為本報告旨在討論多變數多項式求解，如需要隱函數曲面展示，須藉由其他專業軟體方能實現，下面是x^2+y^2+z^2-1=0利用筆者所編之Function 3D Plotter軟體(Marching cubic algorithm)所畫出之3D圖示。
|
上次修改此網站的日期： 2018年12月02日
