[Clinri]

Дана программа решающая нелинейное диференциальное уравнение 2-го порядка, методом Рунге-Кутта, записывающая координаты в массив

Код:

[Выделить все]

Sub RK() 'решение НЛДУ методом Рунге-Кутта

ThisDrawing.SendCommand "_.erase" & vbCr & "_all" & vbCr & vbCr 'очистка
Dim H As Double
Dim hShva As Double
Dim V As Double
Dim W As Double
Dim A As Double
Dim B As Double
Dim F As Double
Dim Z As Double
Dim y As Double
Dim i As Integer
Dim j As Integer

x0 = 0: x = x0 'X0
hShva = 1.5 'высота усиления
W = hShva: y = W  'начальное уловие Y(X0)=Y0
V = 0: Z = V 'начальное уловие dY/dX0
Const N = 50 'число шагов
B = 5 'ширина шва
H = (B - x0) / N 'шаг

Dim plineObj As AcadLWPolyline
Dim points(0 To N * 2 - 1) As Double
j = 0
For i = 1 To N
    F = ZZZ(y, Z): A = H * F: x = x + H / 2
    y = W + V * H / 2 + A * H / 8: Z = V + A / 2
    F = ZZZ(y, Z): B = H * F: Z = V + B / 2
    F = ZZZ(y, Z): C = H * F: x = x + B / 2
    y = W + H * V + H * C / 2: Z = V + C: F = ZZZ(y, Z)
    y = W + H * (V + (A + B + C) / 6): W = y
    Z = V + (A + (B + C) * 2 + H * F) / 6: V = Z
    
    If x < 0 Then Exit For
    If y < 0 Then Exit For
    If Abs(x) > 50 Then Exit For
    points(j) = x: points(j + 1) = y
    j = j + 2
    Debug.Print "для x=" + Str(x) + " y=" + Str(y) + " Dy/Dx=" + Str(Z)
Next i
Debug.Print i

Set plineObj = ThisDrawing.ModelSpace.AddLightWeightPolyline(points)


Dim Center(0 To 2) As Double
Center(0) = 0
Center(1) = 2
Center(2) = 0
ZoomCenter Center, 8
End Sub

Function ZZZ(Y1 As Double, Z1 As Double) As Double
Dim F1 As Double
ak = 4
R0 = 5
F1 = (Y1 / ak ^ 2 - 1 / R0) * (1 + Z1 ^ 2) ^ (3 / 2)
ZZZ = F1
End Function

Отладочное окно:

Код:

[Выделить все]

для x= 4.46868600424883E-02 y= 1.49946870733436 Dy/Dx=-1.06267068050847E-02
для x= 8.93685973602525E-02 y= 1.49787431709848 Dy/Dx=-2.12636613776988E-02
для x= .134040071700682 y= 1.49521529081201 Dy/Dx=-3.19211465760304E-02
для x= .178696107582384 y= 1.49148905848088 Dy/Dx=-.042609515764925
для x= .223331476272288 y= 1.48669200974986 Dy/Dx=-5.33392288879675E-02
для x= .267940877264601 y= 1.48081948140103 Dy/Dx=-6.41208895372246E-02
для x= .312518919082744 y= 1.47386574108045 Dy/Dx=-7.49652833788884E-02
для x= .357060099214 y= 1.46582396709879 Dy/Dx=-8.58834183159451E-02
для x= .401558782971757 y= 1.45668622411181 Dy/Dx=-9.68865667987214E-02
для x= .446009181061615 y= 1.44644343444394 Dy/Dx=-.107986310731544
для x= .490405325604597 y= 1.43508534477097 Dy/Dx=-.119194589469825
для x= .534741044342772 y= 1.42260048782587 Dy/Dx=-.130523751457977
для x= .579009932718809 y= 1.40897613873332 Dy/Dx=-.141986610126296
для x= .623205323480348 y= 1.39419826551287 Dy/Dx=-.153596504746383
для x= .667320253411448 y= 1.37825147321582 Dy/Dx=-.165367367042303
для x= .711347426734979 y= 1.36111894107571 Dy/Dx=-.177313794471693
для x= .755279174659999 y= 1.34278235195389 Dy/Dx=-.189451131231222
для x= .799107410464208 y= 1.32322181324819 Dy/Dx=-.201795558209149
для x= .84282357940075 y= 1.30241576830054 Dy/Dx=-.214364193310085
для x= .88641860259701 y= 1.2803408971851 Dy/Dx=-.227175203821155
для x= .929882813965973 y= 1.2569720055773 Dy/Dx=-.240247932783972
для x= .973205888972222 y= 1.23228190019026 Dy/Dx=-.253603041695205
для x= 1.01637676387712 y= 1.20624124901182 Dy/Dx=-.267262672295305
для x= 1.05938354382158 y= 1.17881842427359 Dy/Dx=-.281250630739524
для x= 1.10221339777756 y= 1.14997932572236 Dy/Dx=-.295592598102765
для x= 1.14485243799516 y= 1.11968718132959 Dy/Dx=-.310316371982138
для x= 1.1872855810701 y= 1.08790232205022 Dy/Dx=-.325452144970367
для x= 1.22949638712913 y= 1.05458192660355 Dy/Dx=-.341032827034041
для x= 1.27146687284299 y= 1.01967973146963 Dy/Dx=-.357094420415436
для x= 1.31317729298031 y= .983145700336511 Dy/Dx=-.373676457680962
для x= 1.35460588394759 y= .944925646048783 Dy/Dx=-.390822516091542
для x= 1.39572856113465 y= .904960796633821 Dy/Dx=-.408580824743382
для x= 1.43651855978512 y= .86318729513476 Dy/Dx=-.42700498515736
для x= 1.47694600637712 y= .819535620648371 Dy/Dx=-.446154831505494
для x= 1.51697740390835 y= .773929915001839 Dy/Dx=-.466097463903526
для x= 1.55657500972016 y= .726287195702419 Dy/Dx=-.486908497800876
для x= 1.59569607812506 y= .676516430879686 Dy/Dx=-.508673585359601
для x= 1.63429193148575 y= .624517445525328 Dy/Dx=-.531490282121227
для x= 1.67230681160735 y= .570179619877559 Dy/Dx=-.555470356091613
для x= 1.70967644698394 y= .513380329524468 Dy/Dx=-.580742669403067
для x= 1.74632624854216 y= .45398306159936 Dy/Dx=-.607456809101858
для x= 1.78216901393511 y= .391835120680809 Dy/Dx=-.63578770970784
для x= 1.81710197333116 y= .326764809263548 Dy/Dx=-.665941605866591
для x= 1.85100294037915 y= .258577927255743 Dy/Dx=-.698163794295922
для x= 1.88372522827909 y= .187053377201425 Dy/Dx=-.732748895631736
для x= 1.91509083224416 y= .111937577854484 Dy/Dx=-.770054630682764
для x= 1.94488113147848 y= 3.29372638447495E-02 Dy/Dx=-.810520633487099
 48

Результат:


собственно вопрос: почему в самом конце точка уходит в координаты (0;0) если в отладочном окошке при этом такого значения не записывалось, или в массиве незаписанные элементы имеют значение 0 ?

[gomer]

Код:

[Выделить все]

If x < 0 Then MsgBox "oops x"
    If y < 0 Then MsgBox "oops y"
    If x < 0 Then Exit For
    If y < 0 Then Exit For

надеюсь объяснять не надо?

[Clinri]

Сообщение от gomer Код: [Выделить все] If x < 0 Then MsgBox "oops x" If y < 0 Then MsgBox "oops y" If x < 0 Then Exit For If y < 0 Then Exit For надеюсь объяснять не надо?

выход из цикла при х=0 у=0 происходит для того чтобы график в первом квадранте рисовался
не понял при чём тут MsgBox "oops y"

хотя я подозреваю что размер массива с запасом взят, и не записанные элементы являются нулями, из за чего и рисуется последняя точка 0;0
хотя не уверен, пока не получилось... дорабатываю с применением ReDim Precerve чтобы размер массива точно соотвествал числу записанных элементов

[gomer]

Сообщение от Clinri хотя я подозреваю что размер массива с запасом взят, и не записанные элементы являются нулями, из за чего и рисуется последняя точка 0;0

[Clinri]

Код:

[Выделить все]

Sub RK() 'решение НЛДУ методом Рунге-Кутта

ThisDrawing.SendCommand "_.erase" & vbCr & "_all" & vbCr & vbCr 'очистка
Dim H As Double
Dim N As Integer
Dim hShva As Double
Dim V As Double
Dim W As Double
Dim A As Double
Dim B As Double
Dim F As Double
Dim Z As Double
Dim y As Double
Dim i As Integer
Dim j As Integer

x0 = 0: x = x0 'X0
hShva = 1.5 'высота усиления
W = hShva: y = W  'начальное уловие Y(X0)=Y0
V = 0: Z = V 'начальное уловие dY/dX0
N = 50 'число шагов
B = 5 'ширина шва
H = (B - x0) / N 'шаг

Dim plineObj As AcadLWPolyline
Dim p() As Double
ReDim p((N * 2) - 1)
j = 0
For i = 1 To N
    F = ZZZ(y, Z): A = H * F: x = x + H / 2
    y = W + V * H / 2 + A * H / 8: Z = V + A / 2
    F = ZZZ(y, Z): B = H * F: Z = V + B / 2
    F = ZZZ(y, Z): C = H * F: x = x + B / 2
    y = W + H * V + H * C / 2: Z = V + C: F = ZZZ(y, Z)
    y = W + H * (V + (A + B + C) / 6): W = y
    Z = V + (A + (B + C) * 2 + H * F) / 6: V = Z
    
    If x < 0 Then Exit For
    If y < 0 Then Exit For
    If Abs(x) > 50 Then Exit For
    
    'записываем первый элемент со значением (x0;y)
        If j = 0 Then
            p(j) = x0: p(j + 1) = y
            Debug.Print "для x"; j; "="; p(j); " y"; j + 1; "="; p(j + 1); " Dy/Dx= 0"
            j = j + 2
        End If
    p(j) = x: p(j + 1) = y
    Debug.Print "для x"; j; "="; p(j); " y"; j + 1; "="; p(j + 1); " Dy/Dx=" + Str(Z)
    j = j + 2
Next i

'записываем последний элемент со значением (x;0)
p(j) = p(j - 2) + (0 - p(j - 1)) * (p(j - 2) - p(j - 4)) / (p(j - 1) - p(j - 3))
p(j + 1) = 0

ReDim Preserve p(j + 1)

Debug.Print i
Debug.Print j + 1

Set plineObj = ThisDrawing.ModelSpace.AddLightWeightPolyline(p)


Dim Center(0 To 2) As Double
Center(0) = 0
Center(1) = 2
Center(2) = 0
ZoomCenter Center, 8
End Sub

Function ZZZ(Y1 As Double, Z1 As Double) As Double
Dim F1 As Double
ak = 4
R0 = 5
F1 = (Y1 / ak ^ 2 - 1 / R0) * (1 + Z1 ^ 2) ^ (3 / 2)
ZZZ = F1
End Function

точно)) размер массива подогнал и всё отлично! СПАСИБО! ))