Mit der Størmer-Verlet-Methode Differentialgleichungen 2. bzw. n.ter Ordnung lösen

    Diese Seite verwendet Cookies. Durch die Nutzung unserer Seite erklären Sie sich damit einverstanden, dass wir Cookies setzen. Weitere Informationen

    • Mit der Størmer-Verlet-Methode Differentialgleichungen 2. bzw. n.ter Ordnung lösen

      Abt. Mit der Størmer-Verlet-Methode Differentialgleichungen 2. bzw. n.ter Ordnung lösen
      =======================================================
      Das Verfahren, das mit vollem Titel "Newton-Størmer-Verlet leapfrog method" heisst, findet sich in fast allen Physik- und Game-Engines zur Berechnung der Newtonschen Bewegungsgesetze. Mit anderen Worten: Alles was explodiert und dann runterfällt wird mit dieser deutlichen, aber doch nicht zu zeitintensiven Verbesserung des einfachsten aller entsprechenden Algorithmen (Explizites Euler-Verfahren) berechnet. Das nachstehende Programm erlaubt eine Beurteilung der Realitätsnähe, weil auch die analytisch-exakte Lösung der Aufgabe bekannt ist. Fünfstellige Genauigkeit reicht da wohl.
      Gruss
      P.S.: Die Variablen-Rückübergabe an Parameterlisten von Fortran-90-Subroutines ist echt ein Graus!!

      Brainfuck-Quellcode

      1. WindowTitle "Anfangswertproblem einer Differentalgleichung 2. Ordnung lösen nach Størmer-Verlet"
      2. ' ---------------------------------------------------------------------------------------------
      3. ' Aus Fortran90 nach XProfan-11 übersetzt 2014-08 von P. Specht, Wien (Österreich)
      4. ' Keine wie auch immer geartete Gewähr! Verwendung auf alleinige Gefahr de(s|r) Anwender(s|in)!
      5. ' Quelle: http://jean-pierre.moreau.pagesperso-orange.fr/Cplus/stormer.txt
      6. '**********************************************************************************************
      7. '* Differential equation y"=f(x,y,y') by Størmer's method *
      8. ' (erweiterbar auf allg. Verlet type, z.B. y'''[...]=f(x,y,y',y''[,y''',...])
      9. '* --------------------------------------------------------- *
      10. '* SAMPLE RUN (BEISPIEL) *
      11. '* Löse y" = 8*y*y / (1+2*x) from x=0 to x=1 *
      12. proc G :parameters x!,y!,z! '= x, y, y'
      13. return 8*sqr(y!)/(1+2*x!)
      14. endproc
      15. proc F :parameters x!,y!,z! '= x, y, y'
      16. return z!
      17. endproc
      18. '* with initial conditions: x(0)=0, y(0)=1 and y'(0)=-2 *
      19. Declare x![4],y![4],z![4],h! ' Initial conditions:
      20. x![1]= 0
      21. y![1]= 1
      22. z![1]=-2 ' dy/dt
      23. h!=0.01 ' Step
      24. '* Compare to exact solution y = 1/(1+2*x) *
      25. Proc Fx :parameters x! ' Exact Solution:
      26. return 1/(1+2*x!)
      27. endproc
      28. '* --------------------------------------------------------- *
      29. '* Differential equation y"=f(x,y,y') by Stormer's method *
      30. '* --------------------------------------------------------- *
      31. '* X Y Y exact Error *
      32. '* 0.000 1.000000 1.000000 0.0000000000 *
      33. '* 0.010 0.980392 0.980392 0.0000000001 *
      34. '* 0.020 0.961538 0.961538 0.0000000295 *
      35. '* 0.030 0.943396 0.943396 0.0000000457 *
      36. '* 0.040 0.925926 0.925926 0.0000000974 *
      37. '* 0.050 0.909091 0.909091 0.0000001285 *
      38. '* ... ... ... ... *
      39. '* 0.950 0.344866 0.344828 0.0000381695 *
      40. '* 0.960 0.342505 0.342466 0.0000388874 *
      41. '* 0.970 0.340176 0.340136 0.0000396196 *
      42. '* 0.980 0.337878 0.337838 0.0000403406 *
      43. '* 0.990 0.335612 0.335570 0.0000410721 *
      44. '* 1.000 0.333375 0.333333 0.0000418231 *
      45. '* *
      46. '* F90 Release By J-P Moreau, Paris (www.jpmoreau.fr). *
      47. '*************************************************************
      48. WindowStyle 24:font 2:set("decimals",16):Window 0,0-%maxx,%maxy-40
      49. Declare c![4],a1!,a2!,a3!,yex!,k&,er! , x1!,y1!,z1!
      50. a1!= 1.083333333333333333 : a2!=-2*(a1!-1) : a3!= a1!-1
      51. clearclip
      52. print "\n------------------------------------------------------------------------------"
      53. print " Differentialgleichung y''=f(x,y,y') mit Störmer's Methode lösen "
      54. print "------------------------------------------------------------------------------\n"
      55. print " X Y Y exact Error\n"
      56. putclip " X | Y | Y_exact | Error \n"
      57. yex!=Fx(x![1]):print tab(2);x![1],tab(22);y![1],tab(42);yex!,tab(62);er!
      58. putclip str$(x![1])+"|"+str$(y![1])+"|"+str$(yex!)+"|"+str$(er!)+"\n"
      59. ' Runge-Kutta for first 2 steps:
      60. whileloop 1,2:k&=&Loop
      61. RK4D2 x![k&],y![k&],z![k&],h!
      62. x![k&+1]=x1!
      63. y![k&+1]=y1!
      64. z![k&+1]=z1!
      65. yex!=Fx(x![k&+1])
      66. er!=abs(yex!-y![k&+1])
      67. print tab(2);x![k&+1],tab(22);y![k&+1],tab(42);yex!,tab(62);er!
      68. putclip str$(x![k&+1])+"|"+str$(y![k&+1])+"|"+str$(yex!)+"|"+str$(er!)+"\n"
      69. endwhile
      70. REPEAT 'Main loop G10:
      71. whileloop 2,4:k&=&Loop
      72. c![k&]=G( x![5-k&],y![5-k&],z![5-k&] )
      73. endwhile
      74. y![4]=2*y![3]-y![2]+sqr(h!)*(a1!*c![2]+a2!*c![3]+a3!*c![4])
      75. x![4]=x![3]+h!
      76. yex!=Fx(x![4])
      77. er!=abs(yex!-y![4])
      78. print tab(2);x![4],tab(22);y![4],tab(42);yex!,tab(62);er!
      79. putclip str$(x![4])+"|"+str$(y![4])+"|"+str$(yex!)+"|"+str$(er!)+"\n"
      80. Whileloop 1,3:k&=&Loop 'für nächsten Schritt umkopieren:
      81. x![k&]=x![k&+1]
      82. y![k&]=y![k&+1]
      83. endwhile
      84. if %csrlin>40:waitinput 4000:cls
      85. print "\n X Y Y-exakt Error\n"
      86. endif
      87. UNTIL x![3]>1 'end x value = 1
      88. print "\n EOF"
      89. Print " Resultate in Zwischenablage!"
      90. Waitinput
      91. END
      92. PROC RK4D2 :parameters x!,y!,z!,h! ',x1!,y1!,z1!
      93. declare c1!,d1!,c2!,d2!,c3!,d3!,c4!,d4!
      94. c1!=F(x!,y!,z!)
      95. d1!=G(x!,y!,z!)
      96. c2!=F(x!+h!/2, y!+h!/2*c1!, z!+h!/2*d1!)
      97. d2!=G(x!+h!/2, y!+h!/2*c1!, z!+h!/2*d1!)
      98. c3!=F(x!+h!/2, y!+h!/2*c2!, z!+h!/2*d2!)
      99. d3!=G(x!+h!/2, y!+h!/2*c2!, z!+h!/2*d2!)
      100. c4!=F(x!+h!, y!+h!*c3!, z!+h!*d3! )
      101. d4!=G(x!+h!, y!+h!*c3!, z!+h!*d3! )
      102. x1!=x!+h!
      103. y1!=y!+h!*(c1!+2*c2!+2*c3!+c4!)/6
      104. z1!=z!+h!*(d1!+2*d2!+2*d3!+d4!)/6
      105. endproc
      106. PROGEND
      107. '{=============================================================================================
      108. ' Stormer, aka Newton-Størmer-Verlet leapfrog method
      109. ' =============================================================================================
      110. ' Namensgebung: Der Algorithmus selbst geht zurück auf Delambre (1791), oftmals wiederentdeckt,
      111. ' u.a. 1907 durch Carl Størmer, 1909 durch Cowell und Cromlin, 1967 populär durch Loup Verlet.
      112. ' Verbessert die explizite Euler-Methode ohne allzu hohen Mehraufwand. Oft in Physik-Engines
      113. ' zu finden. Details siehe http://en.wikipedia.org/wiki/Verlet_integration
      114. ' =============================================================================================
      115. ' Methode STORMER (aka Newton-Störmer-Verlet-leapfrog method):
      116. ' löst Y"=f(x,y,y') mit gegebenen Anfangsbedingungen
      117. ' ----------------------------------------------------------------------------------
      118. ' Eine Differentialgleichung 2. Ordnung (Y") kann ersetzt werden durch ein
      119. ' Differentialgleichungs-SYSTEM bestehend aus 2 Gleichungen 1. Ordnung!
      120. ' Wenn also y"=f(x,y,y') gegeben ist mit y(a=Anfang) und y'(a=Anfang), dann kann man
      121. ' diese Formel durch Substitution u = y' umformen zu:
      122. '
      123. ' | u'=f(x,y,u)
      124. ' | y'=u
      125. ' | mit den neuen Anfangsbedingungen y(a), u(a)
      126. '
      127. ' Beginnen wir die Betrachtung mit jener speziellen Form der Taylor-Entwicklung,
      128. ' bei der der Rest als Integral dargestellt wird (Wiki: Taylorentwicklung)
      129. ' x+h
      130. ' y(x+h) = y(x) + h * y'(x) + INTEGRAL (x+h-t)*y"(t).dt
      131. ' x
      132. ' Durch Verbringung der Nichtintegral-Terme auf die linke Seite und Aufspaltung
      133. ' des 'Restfehler-Integrals' in zwei Teile kann diese Form auch folgender-
      134. ' maßen dargestellt werden:
      135. ' x+h x-h
      136. ' y(x+h) - 2y(x) + y(x-h) = INTEGRAL ... + INTEGRAL ...
      137. ' x x
      138. ' Für das zweite Integral wird die Variable x nun substituiert durch: u = 2*x-t
      139. ' Dann gilt auf Grund der Kettenregel:
      140. ' x-h x+h
      141. ' INTEGRAL (x-h-t) * y"(t).dt = - INTEGRAL (x+h-u)*y"(u).du
      142. ' x x
      143. ' und die Ableitung .du rück-substituiert (zu .-dt) liefert dann:
      144. ' x+h
      145. ' du=-dt ==> = INTEGRAL (x+h-t)*y"(2x-t).dt
      146. ' x
      147. ' Somit kann man schreiben: x+h
      148. ' y(x+h) - 2*y(x) + y(x-h) = INTEGRAL (x+h-t)*[y"(t)+y"(2x-t)].dt
      149. ' x
      150. ' Nun verwenden wir unter Hinweis, daß ja y"(x)=f(x,y,y') gilt, folgenden Ausdruck
      151. ' als Interpolationspolynom:
      152. ' x_n+1
      153. ' y_n+1 - y_n + y_n-1 = INTEGRAL (x_n+1 - t)*[P(t)+P(2*x_n - 1)].dt = ...
      154. ' x_n
      155. ' x_n+1
      156. ' ...= h^2*INTEGRAL (x_n+1 -t)*[O_00*Div(f_n)+O_11*Div(f_n)+O_22*Div(f_n)].dt ,
      157. ' x_n
      158. ' 1 |(-s) (m)|
      159. ' wobei O_mm = -1^m * INTEGRAL (1-s)*| + | .ds
      160. ' 0 |( m) (s)|
      161. ' (vgl. auch 'Adams-Bashford-Methode')
      162. ' Letzteres bringt uns schließlich zu STORMER's Formeln:
      163. ' ==================================================================================
      164. ' Explizit formuliert: y_n+1 -2*y_n + y_n-1 = h^2/12 * [13 * f_n - 2*f_n-1 + f_n-2]
      165. ' Das interessierende y_n+1 = h^2/12 * [13 * f_n - 2*f_n-1 + f_n-2]+ 2*y_n - y_n-1
      166. ' (f_n+1 muss nicht bekannt sein: Prognoserechnung bei aktuell ergänzten Zeitreihen)
      167. ' ----------------------------------------------------------------------------------
      168. ' oder implizit: y_n+1 -2*y_n + y_n-1 = h^2/12 * [f_n+1 + 10 * f_n + f_n-1 ]
      169. '}==================================================================================
      Alles anzeigen
      Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3