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

    Unsere Datenschutzerklärung wurde aktualisiert. Mit der Nutzung unseres Forums akzeptierst Du unsere Datenschutzerklärung. Du bestätigst zudem, dass Du mindestens 16 Jahre alt bist.

    • 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