Differerentialgleichung mittels Prädiktor-Korrektor-Methode 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.

    • Differerentialgleichung mittels Prädiktor-Korrektor-Methode lösen

      Abt. Differerentialgleichung mittels Prädiktor-Korrektor-Methode lösen
      ===========================================
      Ist die Differentialgleichung lediglich von 1. Ordnung, determinieren also z.B. nur aktuelle Augenblickswerte die Änderungsgeschwindigkeit eines Systems (Beispiel: Radioaktiver Zerfall), dann bietet das Prädiktor-Korrektor-Verfahren eine rasch konvergierende Lösung zur schrittweisen Berechnung der Systemverhaltens. Bestimmte Anfangsbedingungen müssen dann allerdings wieder vorgegeben werden. Man spricht vom "Anfangswertproblem", das deutlich leichter lösbar ist als die Frage, welche Funktion stets bestimmte Randwerte/Randbedinungen erfüllt (Randwertproblem, kniffeliger und oft besser mit analytischen Überlegungen lösbar - manchmal aber auch garnicht).
      Gruss

      P.S.: Nachstehend wieder eine Rückübersetzung aus Fortran-90 bzw. ursprünglich Turbopascal. F90 bedient sich bei Subroutinen des Call by Reference, sodaß Parameter-Arithmetik sich auf das Aufrufende Programm auswirkt: Große Stolperfalle!!!

      Brainfuck-Quellcode

      1. WindowTitle " Anfangswert-bestimmte Differentialgleichung "+\
      2. "1. Ordnung mit PRÄDIKTOR-KORREKTOR-METHODE lösen"
      3. '{ (T) Translated to XProfan-11 2014-08 by P.Specht, Wien
      4. '***********************************************************
      5. '* SOLVING DIFFERENTIAL EQUATIONS OF ORDER 1 *
      6. '* (Prediction-correction method) *
      7. '* ------------------------------------------------------- *
      8. ' Beisp. für eine DGL 1. Ordnung: Etwa das Zerfallsgesetz *
      9. ' d.Radioaktivität. Nachstehend: Willkürliche Testfunktion *
      10. '* ------------------------------------------------------- *
      11. '* Reference: *
      12. '* "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc *
      13. '* DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991" *
      14. '* F90 version by J-P Moreau, Paris (www.jpmoreau.fr) *
      15. '* ------------------------------------------------------- *
      16. '* SAMPLE RUN: *
      17. '* (Integrate y'=(4xy+4x*sqrt(y))/(1+x²) from x=0 to x=10 *
      18. '* with initial condition y(0)=1) *
      19. '* *
      20. '* Input x begin: 0 *
      21. '* Input x end : 10 *
      22. '* Input y begin: 1 *
      23. '* Input number of points: 10 *
      24. '* Input finesse: 10 *
      25. '* *
      26. '* X Y *
      27. '* ----------------------------- *
      28. '* 1.000000 8.999984 *
      29. '* 2.000000 80.999881 *
      30. '* 3.000000 360.999496 *
      31. '* 4.000000 1088.998512 *
      32. '* 5.000000 2600.996484 *
      33. '* 6.000000 5328.992838 *
      34. '* 7.000000 9800.986875 *
      35. '* 8.000000 16640.977767 *
      36. '* 9.000000 26568.964559 *
      37. '* 10.000000 40400.946171 *
      38. '* ----------------------------- *
      39. '* *
      40. '}**********************************************************
      41. '{ PROGRAM EQUDIFPC
      42. WindowStyle 24:Font 2:randomize:Window 0,0-%maxx,%maxy-40
      43. declare tx![100],ty![100],x!,xi!,xf!,yi!, fi&,n&,i&
      44. if 1 ' 0=Manual input, 1=Example
      45. xi!=0 :print "\n\n X_begin: ";xi!
      46. xf!=10 :print " X_end : ";xf!
      47. yi!=1 :Print " Y_begin: ";yi!
      48. n& =10 :print " Number of points inbetween: ";n&
      49. fi&=10 :print " Granularity between points: ";fi&
      50. else
      51. print "\n\n Input x begin: "; :input xi!
      52. print " Input x end : "; :input xf!
      53. Print " Input y begin: "; :input yi!
      54. print " Input number of points: ";:input n&
      55. print " Input Granularity: "; :input fi&
      56. endif
      57. equadiff_pc tx![],ty![],xi!,xf!,yi!,n&,fi&
      58. print ' Results:
      59. print " X Y "
      60. print " -----------------------------"
      61. WhileLoop n&:i&=&Loop
      62. print format$("%g",tx![i&]),format$("%g",ty![i&])
      63. endwhile
      64. print " ----------------------------- "
      65. waitinput
      66. END
      67. '} End of main program
      68. ' Hier die Definition der Differentialgleichung f:
      69. ' f: y'=( 4xy + 4x*sqrt(y))/(1+x^2)
      70. proc f :parameters x!,y!
      71. declare f!
      72. if y!>=0
      73. f! = ( 4*x!*y!+4*x!*sqrt(y!) ) / ( 1+sqr(x!) )
      74. endif
      75. return f!
      76. endproc
      77. proc runge_kutta :parameters h!,x!,y!
      78. ' classical Runge-Kutta method of order 4
      79. declare a!,b!,c!,d!,f!
      80. a!=h!*f(x!,y!)
      81. b!=h!*f(x!+h!/2, y!+a!/2)
      82. c!=h!*f(x!+h!/2, y!+b!/2)
      83. d!=h!*f(x!+h!,y!+c!)
      84. y!=y!+(a!+2*b!+2*c!+d!)/6
      85. z!=y!
      86. endproc
      87. '**********************************************************
      88. '{ Prediction-correction method *
      89. '* ------------------------------------------------------ *
      90. '* INPUTS: xi begin x value *
      91. '* xf end x value *
      92. '* y1 begin y value (at x=xi) *
      93. '* m number of points to calculate *
      94. '* fi finesse (number of intermediate *
      95. '* points (for example 10) *
      96. '* OUTPUTS: *
      97. '* tx table of x values (m values) *
      98. '* ty table of y values (m values) *
      99. '* *
      100. '* DESCRIPTION: *
      101. '* The algorithm has the following steps: *
      102. '* 1. calculate y2, y3, y4 using a classical Runge- *
      103. '* Kutta method of order 4. *
      104. '* 2. from point (x4, y4), first estimate y(n+1) by *
      105. '* PREDICTION formula: *
      106. '* y(n+1)=y(n)+h/24(55y'(n)-59y'(n-1)+37y'(n-2)-9y'(n-3) *
      107. '* *
      108. '* then continue with CORREKTION formula: *
      109. '* y(n+1)=y(n) + h/24(9y'(n+1)+19y'(n)-5y'(n-1)+y'(n-2) *
      110. '* *
      111. '* substituting the y'(n+1) = f(x(n+1),y(n+1)) by the *
      112. '* estimated value of y(n+1), until convergence *
      113. '* is obtained. *
      114. '}*********************************************************
      115. Proc equadiff_pc :parameters tx![],ty![],xi!,xf!,yi!,m&,fi&
      116. declare ni&,h!,w!,x!,y!,z!,p![3] , k&,i&,j&
      117. z!=yi! : case (m&>100) or (fi&<1) : Goto "retour"
      118. h!=(xf!-xi!)/fi&/m&
      119. p![3]=f(xi!,yi!) : tx![0]=xi! : ty![0]=yi!
      120. k&=0
      121. whileloop m&:i&=&Loop : ni&=(i&-1)*fi&-1
      122. whileloop fi&:j&=&Loop : x!=xi!+h!*(ni&+j&) : inc k&
      123. if k&<4
      124. runge_kutta h!,x!,z!
      125. x!=x!+h! '<<< wird in Fortran90 in runge_kutta selbst als Rückparameter erledigt!
      126. p![3-k&]=f(x!,z!)
      127. else
      128. x!=x!+h!
      129. w!=z! +h!/24*(55*p![0]-59*p![1] +37*p![2] -9*p![3])
      130. ' y(n+1)=y(n)+h/24*(55*y'(n)-59y'(n-1)+37*y'(n-2)-9*y'(n-3)
      131. Repeat
      132. y!=w!
      133. w!=z!+h!/24*(9*f(x!,y!)+19*p![0]-5*p![1] + p![2])
      134. ' y(n+1)=y(n)+h/24(9*y'(n+1) +19*y'(n)-5*y'(n-1)+y'(n-2)
      135. Until abs(y!-w!)<val("1e-10")
      136. z!=w! : p![3]=p![2] : p![2]=p![1] : p![1]=p![0]
      137. p![0]=f(x!,z!)
      138. endif
      139. endwhile 'j loop
      140. tx![i&]=x!
      141. ty![i&]=z!
      142. endwhile 'i loop
      143. endproc
      Alles anzeigen
      Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3