Algorithmen Teil VIII: Ein Bug kommt selten allein...

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.

  • Algorithmen Teil VIII: Ein Bug kommt selten allein...

    Und wieder ein neues Kapitel verspielten Unsinns, fehlgeschlagener Versuche, haarstreubender Ergebnisse und unwiederholbarer Fehler. Viel Theorie und gar wenig Praxis-taugliches wartet wieder auf möglichst kreatives Recycling. Willkommen also, ihr edlen BesucherInnen unserer (Software-)Mülldeponie, - möge der Stockschnupfen mit Euch sein!

    Abt. Mut du ma Ordnung!
    ================
    Was waren bisher eigentlich die Themen? Kann man das mal irgendwie in eine Ordnung bringen? Ich versuch's mal.
    Gruss
    Empirisch ermittelte Themenblöcke aus den bisherigen Kapiteln (=Threads) I.-VII

    A. INFORMATIK ALLG., inkl. Methoden & Paradigmen

    B. ALGORITHMEN
    B.1 Stringalgorithmen
    B.2 Sortieren
    B.3 Symbolverarbeitung und KI
    B.4 Graphikroutinen
    B.5 Optimierung in Technik und Wirtschaft

    C. DATENHALTUNG, DATENBANKEN
    D. DATENKOMPRESSION
    E. DATENSICHERHEIT

    F. XPROFAN: EIGENHEITEN, TIPPS & TRICKS

    G. WINDOWS
    G.1 Windows-API

    H. INTERNET
    I. UTILITIES
    J. BENCHMARKING
    K. ASSEMBLER
    L. HARDWARE

    M. MATHEMATIK & LOGIK
    M.1 ZAHLENTHEORIE
    M.2 KOMBINATORIK
    M.3 ALGEBRA UND NUMERIK
    M.4 MATRIZENRECHNUNG
    M.5 DIFFERENTIALGLEICHUNGSSYSTEME
    M.6 GEOMETRIE
    M.9 UNTERHALTUNGSMATHEMATIK

    N.1 PHYSIK
    N.2 ASTRONOMIE, RAUMFAHRT

    O. TECHNIK & MEDIZIN

    P. MULTIMEDIA
    P.1 VIDEO
    P.2 AUDIO
    P.3 SPIELE

    Q. WIRTSCHAFT

    R. KUNST, KULTUR, GESELLSCHAFT
    S. SONSTIGES: PHILOSOPHISCHES, KURIOSES, HUMOR

    --------------
    P.S.: ... und fast überall war Unsinn drin!
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3
  • Abt. Spezielle Datenhaltung: Die DEQue-Struktur
    ===============================
    Erläuterungen im Programmtext!
    Gruss

    P.S.: Zit.:'In der Praxis verwendet man die Deque unter anderem ... zur Textsuche mittels regulärer Ausdrücke (Pattern-Matching-Algorithmus)'

    Quellcode

    1. WindowTitle "Datenstruktur DEQue (DoubleEnded Queue): Die beidseitige Warteschlange (aka Doppelstack)"
    2. '(CL) CopyLeft 2014-08 by P. Specht, Vienna/Austria
    3. ' Keine wie auch immer geartete Gewähr, Verwendung auf alleinige Gefahr der AnwenderInnen!
    4. WindowStyle 24:CLS 'Window 0,0-%maxx,%maxy-40
    5. Font 2
    6. Print "\n============================================================"
    7. print " KLEINES DEQue-TUTORIAL "
    8. print " PROCs are Freeware, Tutorial=CopyLeft by P. Specht "
    9. Print "============================================================"
    10. Print "\n Beim Umsetzen z.B. des Melkman Convex Hull-Algorithmus [Melkman 1987] "
    11. Print " begegnet man einer Datenstruktur namens DEQue (spr.'Deck' Double Ended Queue)"
    12. print " Während 'Queue' im Englischen 'Warteschlange' bedeutet (aka. FIFO-Speicher)"
    13. print " handelt es sich beim DEQue um eine Warteschlange, bei der zB. Privatpatien-"
    14. print " ten vorgezogen werden können, aber auch Leute die frustrierend lange warten"
    15. print " sich hinten wieder entfernen (oder die Schlange wechseln) können. "
    16. print "\n Auf Datenhaltungsstrukturen bezogen bedeutet die Realisierung, daß entweder"
    17. print " Arrayelemente physisch stellenweise vor- bzw. zurückverschoben werden (hier"
    18. print " im Programm realisiert, relativ langsam), oder vor-&rück-verlinkte Listen "
    19. print " die entsprechende Funktionalität bereitstellen (Aufwendig, am besten als "
    20. print " Memory-Verwaltungssystem auf maschinennaher Ebene, Betriebssystemabhängig)."
    21. print "\n Hinweis: Einfache Warteschlangen benötigten nur die Hälfte der PROCs hier! "
    22. print "\n Der Begriff TOP bedeutet hier 'OBEN', der Begriff Btm=Bottom ='UNTEN' "
    23. print "\n Quelle: http://de.wikipedia.org/wiki/Deque [Taste]"
    24. Font 0
    25. waitinput
    26. Cls
    27. Declare Wert&,DEQ&[],numw&
    28. numw&=4 'temporary numwidth, mit der das DEQue zu Testzwecken ausgegeben wird
    29. print "\n MakeConstDEQue DEQ&[],11,-7"
    30. MakeConstDEQue DEQ&[],11,-7
    31. ShowDEQue DEQ&[],numw&
    32. Waitinput
    33. print "\n MakeConstDEQue DEQ&[],11,88"
    34. MakeConstDEQue DEQ&[],11,88
    35. ShowDEQue DEQ&[],numw&
    36. Waitinput
    37. print "\n MakeConstDEQue DEQ&[],11 <<< Ohne Wertangabe wird NICHT überschrieben!"
    38. MakeConstDEQue DEQ&[],11
    39. ShowDEQue DEQ&[],numw&
    40. Waitinput
    41. print "\n MakeIndexDEQue DEQ&[],7"
    42. MakeIndexDEQue DEQ&[],7
    43. ShowDEQue DEQ&[],numw&
    44. Waitinput
    45. print "\n MakeRevertIndexDEQue DEQ&[],11"
    46. MakeRevertIndexDEQue DEQ&[],11
    47. ShowDEQue DEQ&[],numw&
    48. Waitinput
    49. skip:
    50. CLS
    51. print "\n Clear DEQ&[]"
    52. Clear DEQ&[]
    53. ShowDEQue DEQ&[],numw&
    54. waitinput
    55. Print "\n Mehrfaches PutTop(Wert&,DEQ&[])"
    56. DEQ&[]=PutTop(555,DEQ&[])
    57. DEQ&[]=PutTop(111,DEQ&[])
    58. DEQ&[]=PutTop(333,DEQ&[])
    59. DEQ&[]=PutTop(777,DEQ&[])
    60. ShowDEQue DEQ&[],numw&
    61. Waitinput
    62. Print "\n Mehrfaches PopTop(DEQ&[]) <<< mehr als Werte in der Schlange stehen)"
    63. font 2:print " Top=";PopTop(DEQ&[]):font 0
    64. ShowDEQue DEQ&[],numw&
    65. Waitinput
    66. font 2:print " Top=";PopTop(DEQ&[]):font 0
    67. ShowDEQue DEQ&[],numw&
    68. Waitinput
    69. font 2:print " Top=";PopTop(DEQ&[]):font 0
    70. ShowDEQue DEQ&[],numw&
    71. Waitinput
    72. font 2:print " Top=";PopTop(DEQ&[]):font 0
    73. ShowDEQue DEQ&[],numw&
    74. waitinput
    75. print " Magic Error-Wert =";int(999999999);" (= minus 9 mal die 9)"
    76. font 2:print " Top=";PopTop(DEQ&[]):font 0
    77. ShowDEQue DEQ&[],numw&
    78. waitinput
    79. font 2:print " Top=";PopTop(DEQ&[]):font 0
    80. ShowDEQue DEQ&[],numw&
    81. waitinput
    82. CLS
    83. Print "\n Clear DEQ&[]"
    84. Clear DEQ&[]
    85. ShowDEQue DEQ&[],numw&
    86. print
    87. waitinput
    88. print "\n Mehrfaches PutBtm(Wert&,DEQ&[])"
    89. DEQ&[]=PutBtm(666,DEQ&[])
    90. ShowDEQue DEQ&[],numw&
    91. waitinput
    92. DEQ&[]=PutBtm(444,DEQ&[])
    93. ShowDEQue DEQ&[],numw&
    94. waitinput
    95. DEQ&[]=PutBtm(222,DEQ&[])
    96. ShowDEQue DEQ&[],numw&
    97. waitinput
    98. print "\n Mehrfaches Wert&=PopBtm(DEQ&[])"
    99. font 2:print " Btm=";PopBtm(DEQ&[]):font 0
    100. ShowDEQue DEQ&[],numw&
    101. waitinput
    102. font 2:print " Btm=";PopBtm(DEQ&[]):font 0
    103. ShowDEQue DEQ&[],numw&
    104. waitinput
    105. font 2:print " Btm=";PopBtm(DEQ&[]):font 0
    106. ShowDEQue DEQ&[],numw&
    107. waitinput
    108. font 2:print " Btm=";PopBtm(DEQ&[]):font 0
    109. ShowDEQue DEQ&[],numw&
    110. waitinput
    111. font 2:print " Btm=";PopBtm(DEQ&[]):font 0
    112. ShowDEQue DEQ&[],numw&
    113. waitinput
    114. CLS
    115. Print "\n Clear DEQ&[]"
    116. CLEAR DEQ&[]
    117. ShowDEQue DEQ&[],numw&
    118. waitinput
    119. Print "\n Mehrfaches PutTop(.) / PutBtm(), wild durcheinander"
    120. DEQ&[]=PutTop(555,DEQ&[])
    121. DEQ&[]=PutBtm(444,DEQ&[])
    122. DEQ&[]=PutTop(666,DEQ&[])
    123. ShowDEQue DEQ&[],numw&
    124. waitinput
    125. Print "\n Mehrfaches PopTop / PopBtm, wild durcheinander"
    126. font 2:print " Top=";PopTop(DEQ&[]):font 0
    127. ShowDEQue DEQ&[],numw&
    128. waitinput
    129. font 2:print " Btm=";PopBtm(DEQ&[]):font 0
    130. ShowDEQue DEQ&[],numw&
    131. waitinput
    132. font 2:print " Btm=";PopBtm(DEQ&[]):font 0
    133. ShowDEQue DEQ&[],numw&
    134. waitinput
    135. font 2:print " Top=";PopTop(DEQ&[]):font 0
    136. ShowDEQue DEQ&[],numw&
    137. waitinput
    138. CLS
    139. Print "\n Clear DEQ&[]"
    140. CLEAR DEQ&[]
    141. ShowDEQue DEQ&[],numw&
    142. waitinput
    143. Print "\n Mehrfaches PutTop(.) / PutBtm(), wild durcheinander"
    144. DEQ&[]=PutTop(555,DEQ&[])
    145. DEQ&[]=PutBtm(444,DEQ&[])
    146. DEQ&[]=PutTop(666,DEQ&[])
    147. ShowDEQue DEQ&[],numw&
    148. waitinput
    149. font 2:print " CopyTop=";CopyTop(DEQ&[]):font 0
    150. ShowDEQue DEQ&[],numw&
    151. waitinput
    152. font 2:print " CopyBtm=";CopyBtm(DEQ&[]):font 0
    153. ShowDEQue DEQ&[],numw&
    154. Print "\n\n============================================================"
    155. print " ENDE DES KLEINEN DEQue-TUTORIALS "
    156. print " Gruss, P. Specht "
    157. Print "============================================================"
    158. waitinput
    159. End
    160. proc MakeIndexDEQue :parameters q&[],n&
    161. SetSize q&[],n&
    162. q&[]=&index+1
    163. return q&[]
    164. endproc
    165. proc MakeRevertIndexDEQue :parameters q&[],n&
    166. SetSize q&[],n&
    167. q&[]=n&-&index
    168. return q&[]
    169. endproc
    170. proc MakeConstDEQue :parameters q&[],n&,wert&
    171. SetSize q&[],n&
    172. q&[]=Wert&
    173. return q&[]
    174. endproc
    175. Proc PopBtm :parameters q&[]
    176. 'pops bottom and returns popped bottom element
    177. var sz&=sizeof(q&[])
    178. if sz&>0
    179. var wert&=q&[0]
    180. q&[]=q&[&index+(&index<(sz&-1))]
    181. setsize q&[],sizeof(q&[])-(sz&>0)
    182. return wert&
    183. else
    184. return int(-999999999)
    185. endif
    186. endproc
    187. proc PopTop :parameters q&[]
    188. declare wert&
    189. ' pops and returns the popped top element
    190. var sz&=sizeof(q&[])
    191. if sz&>0
    192. wert&=q&[sz&-1]
    193. else
    194. wert&=int(-999999999)
    195. endif
    196. setsize q&[],sz&-(sz&>0)
    197. return wert&
    198. endproc
    199. Proc PutTop :parameters Wert&,q&[]
    200. var sz&=sizeof(q&[])
    201. inc sz&
    202. setsize q&[],sz&
    203. q&[sz&-1]=Wert&
    204. endproc
    205. proc PutBtm :parameters Wert&,q&[]
    206. declare tmp&[]
    207. var sz&=sizeof(q&[])
    208. inc sz&
    209. setsize q&[],sz&
    210. tmp&[]=q&[]
    211. q&[]=tmp&[&index-(&index>0)]
    212. clear tmp&[]
    213. q&[0]=Wert&
    214. endproc
    215. Proc CopyBtm :parameters q&[]
    216. ' Returns copy of bottom element
    217. var sz&=sizeof(q&[])
    218. if sz&>0
    219. return q&[0]
    220. else
    221. return int(-999999999)
    222. endif
    223. endproc
    224. proc CopyTop :parameters q&[]
    225. declare wert&
    226. ' returns copy of top element
    227. var sz&=sizeof(q&[])
    228. if sz&>0
    229. wert&=q&[sz&-1]
    230. else
    231. wert&=int(-999999999)
    232. endif
    233. return wert&
    234. endproc
    235. proc ShowDEQue :parameters q&[],nw&
    236. var tmp&=get("numwidth"):set("numwidth",nw&)
    237. var sz&=sizeof(q&[])
    238. print tab(2);"n=";sz&;": ";
    239. whileloop 0,sz&-1
    240. print q&[&Loop],
    241. endwhile
    242. print
    243. set("numwidth",tmp&)
    244. endproc
    Alles anzeigen
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3

    Dieser Beitrag wurde bereits 2 mal editiert, zuletzt von p. specht ()

  • Ergänzende Funktionen zu oben:
    ====================

    Quellcode

    1. proc CopyFromDEQue :parameters q&[],n&
    2. declare wert&
    3. ' returns copy of top element
    4. var sz&=sizeof(q&[])
    5. if (sz&>0) and (n&>=0) and (n&<sz&)
    6. wert&=q&[n&]
    7. else
    8. wert&=int(-999999999)
    9. endif
    10. return wert&
    11. endproc
    12. Proc CmpInDEQue :parameters q&[],n&,m&
    13. 'Compares two elements in DEQue
    14. var sz&=sizeof(q&[])
    15. if (sz&>0) and (n&>=0) and (n&<sz&) and (m&>=0) and (m&<sz&)
    16. return int((q&[n&]>q&[m&])-(q&[n&]<q&[m&]))
    17. else
    18. return int(-999999999)
    19. endif
    20. Endproc
    21. proc SwapInDEQue :parameters q&[],n&,m&
    22. 'Swaps two elements in Queue
    23. var sz&=sizeof(q&[])
    24. if (sz&>0) and (n&>=0) and (n&<sz&) and (m&>=0) and (m&<sz&)
    25. var tmp&=q&[n&]:q&[n&]=q&[m&]:q&[m&]=tmp&
    26. return int(n&<>m&)
    27. else
    28. return int(-999999999)
    29. endif
    30. endproc
    31. proc DelFromDEQue :parameters q&[],n&
    32. 'Delets member n from DEQue
    33. var sz&=sizeof(q&[])
    34. declare wert&
    35. if (sz&>0) and (n&>0) and (n&<sz&)
    36. wert&=q&[n&]
    37. q&[]=q&[&index+(&index>=n&)*(&index<(sz&-1))]
    38. dec sz&
    39. setsize q&[],sz&
    40. else
    41. wert&=int(-999999999)
    42. endif
    43. return wert&
    44. endproc
    Alles anzeigen

    ... mit dem zugehörigen Testteil (Am Ende des Hauptprogramms einfügen):

    Quellcode

    1. Print " Wert& = CopyFromDEQue(DEQ&[],1) "
    2. font 2:print " n.Element=";CopyInDEQue(DEQ&[],1):font 0
    3. ShowDEQue DEQ&[],numw&
    4. waitinput
    5. Print " Wert& = CopyFromDEQue(DEQ&[],11111) "
    6. font 2:print " n.Element=";CopyInDEQue(DEQ&[],11111):font 0
    7. ShowDEQue DEQ&[],numw&
    8. waitinput
    9. CLS
    10. print "\n (n>m)? = CmpInDEQue(q&[],0,1)"
    11. font 2:print " LT-1 EQ0 GT1 =";CmpInDEQue(DEQ&[],0,1):font 0
    12. ShowDEQue DEQ&[],numw&
    13. waitinput
    14. Print "\n (n=m)? = SwapInDEQue(q&[],0,1)"
    15. font 2:print " (N=M)=";SwapInDEQue(DEQ&[],0,1):font 0
    16. ShowDEQue DEQ&[],numw&
    17. waitinput
    18. print "\n (n>m)? = CmpInDEQue(q&[],0,1)"
    19. font 2:print " LT-1 EQ0 GT1 =";CmpInDEQue(DEQ&[],0,1):font 0
    20. ShowDEQue DEQ&[],numw&
    21. waitinput
    22. Print " ok = DelFromDEQue(DEQ&[],1)"
    23. DEQ&[]=DelFromDEQue(DEQ&[],1)
    24. ShowDEQue DEQ&[],numw&
    25. waitinput
    Alles anzeigen
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3

    Dieser Beitrag wurde bereits 1 mal editiert, zuletzt von p. specht ()

  • Ergänzung II zu oben: Ob die gewählten Befehlsbezeichnungen gut zu merken sind, sei mal dahingestellt. Laut Wikipedia gibt es PUT/GET für die Top-Seite und PUSH/POP für die Bottom-Seite, aber die sind bekanntlich schon besetzt. Wer immer damit öfter arbeitet, wird sich vermutlich ein eigenes Schema zulegen. Mein Vorschlag im Programm oben lautete:

    CLEAR q&[] 'natives XProfan
    sz&=SizeOf(q&[]) 'sz=5 bedeutet: 0,1,2,3,4
    SetSize q&[],sizeof(q&[]) 'erweitert das Array also um 1 Stelle
    q&[]=MakeIndexDEQue(q&[],n&)
    q&[]=MakeRevertIndexDEQue(q&[],n&)
    q&[]=MakeConstDEQue(q&[],n&,wert&)
    Wert& = PopBtm(q&[]) '-999999999 = out of range error
    Wert& = PopTop(q&[]) '-999999999 = out of range error
    PutTop Wert&,q&[]
    PutBtm Wert&,q&[]
    Wert& = CopyBtm(q&[])'-999999999 = out of range error
    Wert& = CopyTop(q&[])'-999999999 = out of range error
    ShowDEQue q&[],nw& 'in: Temporary NumWidth
    Wert& = CopyFromDEQue(q&[],n&)'-999999999 = out of range
    LogVgl&=CmpInDEQue(q&[],n&,m&)' n=Pos1,m=Pos2: <:1, '=':0, >: -1
    NisM& = SwapInDEQue(q&[],n&,m&) '(N=M)=1, -999999999=range err
    Wert& = DelFromDEQue(q&[],n&) '-999999999 = out of range

    P.S.: Das RETURN mit Übergabe des Dynamischen Arrays bei den ersten drei Befehlen scheint ausserdem überflüssig zu sein, da im Falle Dynamischer Arrays die Proc-Übergabe per Call by Reference zu erfolgen scheint. Es wird also ohnehin direkt auf das übergebene Array gearbeitet. Bin am klären...
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3
  • Abt. Ist doch logisch, oder?
    =================
    Während sich Programmierer in der Regel in der Welt zweiwertiger Logik (0 - 1, Falsch - Wahr, Aus - Ein) von George Boole und De Morgan bewegen, hält die Realität oft auch verschiedene Zwischenwerte bereit. Man spricht von Mehrwertiger Logik. So haben sich Reichenbach und Lukasiewiczs mit dreiwertiger Logik (z.B. in der Quantenmechanik) und ihren Erweiterungen in Richtung unendlichwertiger Kontinuumslogik (alles zwischen 0 und 1) beschäftigt und dabei deren mathematische und philosophische Konsequenzen erforscht. Kleene und Gödel konnten im Bereich vier- und fünfwertiger Logik Fortschritte erzielen (Kleene untersuchte damit die Grade der Lösbarkeit mathematischer Probleme und erfand dabei die Regular Expressions) und Popper schloss auf dieser Basis auf ethisch vertretbare Mechanismen der verschiedenen wissenschaftlichen Beweisverfahren. Schließlich kehrte das Thema in die Elektronikfertigung und -Qualitätssicherung zurück: In der Hardwarebeschreibungssprache VHDL wird zum Beispiel die im IEEE-Standard 1164 definierte neunwertige Logik verwendet:
    Standard Logic 1164:
    U undefiniert
    X unbekannt (für immer verborgen)
    0 logische Null (sicher falsch)
    1 logische Eins (sicher richtig)
    Z abgekoppelt (hohe Impedanz Z)
    W unbekannt (momentan unbekannt)
    L logische Null low (vermutlich falsch)
    H logische Eins high (möglicherweise richtig)
    - unwichtig, egal
    Man sieht, daß hier die Tendenz bereits Richtung Fuzzy-Logik von L. Zadeh wandert, - mit dem Unterschied, daß dort zusätzlich eine 'innere Interpretation' sehr gutes Sprachgefühl erfordert und der Programmierer diese letztlich erst der Gruppe der Anwender in zeitaufwändigen Gesprächen 'entlocken' muss.
    Gruss

    P.S.: Der guten Ordnung halber sei noch Brouwer's Intuitionistische Logik genannt, die allerdings eher erkenntnistheoretische Grundlagen zum Mechanismus der Forschung legt - Auf Kurzform gebracht: Ohne Bauchgefühl, diffuse Annahmen, wilde Theorien und Widerlegungsversuche zu Vorurteilen gibt's erst mal keinen Erkenntnisgewinn!
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3

    Dieser Beitrag wurde bereits 4 mal editiert, zuletzt von p. specht ()

  • Bugs kommen eben selten allein: KORREKTUR zu vor-voriger Beitrag:
    SetSize q&[],sizeof(q&[])+1 'erweitert das Array also um 1 Stelle!!!
    (SetSize funktioniert anders als Declare )
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3
  • Abt. Wann und wo gehen die Planeten auf?
    ===========================
    Auf der genialen Homepage von Dr. Jean Pierre Moreau findet sich zu dieser Frage ein Programm aus der Zeit vor 1982, dessen Gütligkeit wegen der angewendeten Genauigkeit vermutlich auch 32 Jahre später noch gute Ergebnisse liefert. Ich konnte zwei sinnstörende Fehler im Basic-Originaltext erkennen und beheben, da ich aber leider kein Äquatorial-montiertes Teleskop besitze, kann ich die Ergebnisse derzeit nicht verifizieren. Für entsprechende Kommentare wäre ich natürlich dankbar.
    Gruss

    Brainfuck-Quellcode

    1. WindowTitle upper$(" P o s i t i o n d e r P l a n e t e n b e s t i m m e n ")+" - Achtung, Daten aus 1980!"
    2. ' Source: http://jean-pierre.moreau.pagesperso-orange.fr/Basic/planets_bas.txt
    3. ' Translated to XProfan-11 (CL) CopyLeft 2014-08 by P. Specht, Vienna
    4. ' No warranty whatsoever! Keine wie auch immer geartete Gewähr!
    5. ' Bahndaten aus der Zeit vor 1982, zu verifizieren!
    6. '
    7. '*****************************************************
    8. '* Position der Planeten *
    9. '* ------------------------------------------------- *
    10. '* This program calculates the ephemeris of a planet *
    11. '* in our solar system to adjust an equatorial tele- *
    12. '* scope (Pgm für aequatorial-montierte Teleskope *
    13. '* ------------------------------------------------- *
    14. '* Ref.: "Mathematiques par l'informatique indivi- *
    15. '* duelle - Programmes en BASIC, MASSON, *
    16. '* Paris, 1982, p 100 - 105" . *
    17. '* ------------------------------------------------- *
    18. '* Inputs: *
    19. '* Date: day,month,year *
    20. '* Hour UT: hour *
    21. '* Planet number: 1 to 8 *
    22. '* (Mercury: 1 Venus : 2 Mars : 4 Jupiter: 5 *
    23. '* Saturn : 6 Uranus: 7 Neptune: 8) *
    24. '* Outputs: *
    25. '* Ascent in hours (0 to 24) *
    26. '* Declination in deg. (-90 to 90 North or South) *
    27. '* *
    28. '* SAMPLE RUN: *
    29. '* (Find ascent and declination of planet Mars on *
    30. '* March 10th, 1982 at 6h UT) *
    31. ' *
    32. ' Hinweis von P.Specht:
    33. ' Am 10. März 1982 standen Merkur, Venus, Erde,
    34. ' Mars, Jupiter, Saturn und Pluto alle in einer Linie
    35. ' und auf der selben Seite der Sonne innerhalb eines
    36. ' Winkel von 95 Grad. Die prophezeiten Katastrophen,
    37. ' etwa Erdbeben in Kalifornien sind gsd ausgeblieben.
    38. ' http://www.borderlands.de/net_pdf/NET1107S30-34.pdf
    39. '
    40. '* Date (D,M,Y): 10,3,1982 *
    41. '* Hour UT: 6 *
    42. '* Mercury: 1 Venus : 2 Mars : 4 Jupiter: 5 *
    43. '* Saturn : 6 Uranus: 7 Neptune: 8 *
    44. '* Planet number: 4 *
    45. '* *
    46. '* Ascent : 13 H 8 MN *
    47. '* Declination: 3 ø 45 MN S *
    48. '* *
    49. '* This program allows building the following table: *
    50. '* *
    51. '* Date: March 10th, 1982 at 6H UT. *
    52. '* *
    53. '* Planet Ascent Declination *
    54. '* Mercury 21 H 51 14 ø 45 mn S *
    55. '* Venus 20 H 26 14 ø 57 mn S *
    56. '* Mars 13 H 08 3 ø 45 mn S *
    57. '* Jupiter 14 H 32 13 ø 30 mn S *
    58. '* Saturn 13 H 22 5 ø 42 mn S *
    59. '* *
    60. '* English Version By J-P Moreau, Paris. *
    61. '* (www.jpmoreau.fr) *
    62. '*****************************************************
    63. WindowStyle 24:Window 0,0-%maxx,%maxy:font 2
    64. declare i&,j&,jj&,n&,m&,mm&,aa!,h!,t!,xl!,o!,p!,e!,xi!,xa!,xm!
    65. declare u!,r!,v!,q!,xj!,x!,y!,z!,G!,pi!,xx!,yy!,zz!,d!,b!,c!
    66. declare A![9,8],B![12],A$,data$
    67. 'Init planetary constants (fill table A by columns)
    68. data$="4.01166, 0.071425454, 1.32493, 0.000000742289, 0.823045, "+\
    69. "0.000000566185, 0.205615, 0.122225, 0.387099, "+\
    70. "3.60861, 0.027963119, 2.271616, 0.00000065572, 1.32291,"+\
    71. "0.000000436681, 0.006816, 0.0592301, 0.723332,"+\
    72. "1.72727, 0.0172028, 1.76688, 0.000000818559, 0,"+\
    73. "0, 0.016751, 0, 1,"+\
    74. "2.17756, 0.0091467658, 5.83378, 0.000000879297, 0.851616, "+\
    75. "0.000000371232, 0.093309, 0.0322939, 1.5236,"+\
    76. "4.68279, 0.00145099, 0.2289, 0.000000857, 1.73578,"+\
    77. "0.000000482933, 0.048376, 0.0228418, 5.202799,"+\
    78. "4.8567, 0.00058484, 1.5974, 0.000000412, 1.96856,"+\
    79. "0.000000417308, 0.054311, 0.0435026, 9.552098, "+\
    80. "4.3224, 0.000205424, 2.9523, 0.000000762, 1.2825, "+\
    81. "0.000000238237, 0.047319, 0.013482, 19.21694, "+\
    82. "1.5223, 0.000105061, 0.7637, 0.000000393, 2.28102,"+\
    83. "0.00000052517, 0.008262, 0.0310536, 30.12912"
    84. '=Spaltenweise! Schnell (=in der inneren Schleife) läuft daher der ZEILEN-Index i&!)
    85. whileloop 8:j&=&Loop:whileLoop 9:i&=&loop
    86. A![i&,j&]=val(substr$(data$,(i&+(j&-1)*9),","))
    87. ' print i&;",";j&;"=",int(1+(i&-1)+(j&-1)*9);":",format$("%g",a![i&,j&])," ",:waitinput
    88. endwhile :endwhile ' Basis [1,1] wird hier beibehalten
    89. 'Calendar constants (Ermittlung der Tagesnummer im Jahr)
    90. 'Sinnstörender Komma-Fehler im Originaltext!
    91. 'Schaltjahre sind hier aufgeteilt: Die Erde selbst schaltet ja nicht ;-)
    92. whileloop 12 :B![&Loop]=\
    93. val(substr$("0,31,59.25,90.25,120.25,151.25,181.25,212.25,243.25,273.25,304.25,334.25",&Loop,","))
    94. endwhile
    95. print "\n\n\n "
    96. 'print " Please input date (DD, MM, YYYY) and hour (GMT) of interest:"
    97. print " Bitte das interessierende Datum (DD, MM, YYYY) und die UTC-Stunde (GMT) eingeben:"
    98. print "\n Tag xx:",:input j&
    99. print " Monat xx:",:input mm&
    100. print " Jahr xxxx:",:input aa!
    101. print "\n Stunde UT = GMT(Greenwich): ",:input h!
    102. 'Calculate time t
    103. t!=365.25*(aa!-1901)+B![mm&]+j& 'Tagesnummer seit 1901
    104. t!=INT(t!)+h!/24 'Stunden in Dezimalform dazu
    105. 'calculate earth coordinates
    106. j&=3 : gosub "S1000" 'planet #3 coordinates
    107. g!=x! : h!=y! 'save earth coordinates
    108. Nochmal:
    109. print
    110. print " Merkur/y=1, Venus=2,((Earth=3)), Mars=4, Jupiter=5, Saturn=6, Uranus=7, Neptun=8"
    111. print " Please Input Planet Number: ", :input n&
    112. print
    113. 'calculate coordinates of planet #n
    114. j&=n& : gosub "S1000"
    115. 'calculate geocentric equatorial coordinates
    116. x!=x!-g! : y!=y!-h!
    117. t!=y!*0.917484-z!*0.397772
    118. z!=y!*0.397772+z!*0.917484 : y!=t!
    119. 'calculate ascent and declination
    120. pi!=pi() '=4*arctan(1)
    121. xx!=x! : yy!=y! : gosub "S2000" 'aa=ATN2(y,x) ' Simuliert echte Atn2( , )-Funktion
    122. aa!=zz!
    123. yy!=z! : xx!=SQRT(x!*x!+y!*y!) : gosub "S2000" ' d=ATN2(z, SQR(x*x+y*y) )
    124. d!=zz!
    125. 'conversion
    126. aa!=aa!*12/pi!
    127. case aa!<0 :aa!=24+aa!
    128. h!=INT(aa!) : m&=INT(60*(aa!-h!))
    129. d!=d!*180/pi!
    130. A$="N":case d!<0:A$="S"
    131. d!=ABS(d!) : b!=INT(d!) : c!=INT(60*(d!-b!))
    132. 'print results
    133. print
    134. print "\n Ascension (Aufgangszeit): ";format$("00",h!);" Uhr ";format$("00",m&);" min"
    135. print "\n bei Deklination (Fusspunkt-Winkel): ";
    136. print format$("#00",b!);"° ";format$("00",c!);" ";A$
    137. waitinput
    138. ' Auf Grund der unglücklichen Zweitverwendung von Variablen im Original stets Neustart erforderlich!
    139. END
    140. S1000:
    141. 'planet coordinates
    142. 'calculate planetary constants
    143. XL!=A![1,j&] + A![2,j&] * t!
    144. O! =A![3,j&] + A![4,j&] * t!
    145. P! =A![5,j&] + A![6,j&] * t!
    146. E! =A![7,j&]
    147. XI!=A![8,j&]
    148. XA!=A![9,j&]
    149. 'solve Kepler's equation
    150. xm! = XL!-O! : u!=xm!
    151. whileloop 100:j&=&Loop
    152. u! = xm! + E! * sin(u!)
    153. endwhile
    154. 'formula (3) of reference book
    155. r!=XA!*(cos(u!)-E!)
    156. v!=XA!*SQRT(1-E!*E!)*sin(u!)
    157. O!=O!-P! : xm!=sin(O!) : O!=cos(O!)
    158. q!=sin(P!) : P!=cos(P!)
    159. xj!=sin(XI!) : XI!=cos(XI!)
    160. x!=(P!*O!-XI!*q!*xm!)*r!+(-P!*xm!-XI!*q!*O!)*v!
    161. y!=(q!*O!+XI!*P!*xm!)*r!+(-q!*xm!+XI!*P!*O!)*v!
    162. z!=xj!*xm!*r!+xj!*O!*v!
    163. return
    164. S2000:
    165. 'calculate zz = ATN2(yy,xx)
    166. if ABS(xx!)<val("1e-12")
    167. if yy!>0 : zz!=Pi!/2
    168. elseif ABS(yy!)<val("1e-12"): zz!=-PI!/2
    169. else : zz!=0
    170. endif
    171. else
    172. zz!=Arctan(yy!/xx!)
    173. case xx!<0 :zz!=zz!+Pi!
    174. endif
    175. return
    176. 'end of file planets.prf
    Alles anzeigen
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3
  • 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
  • Anmerkung zu oben: Carl Størmer, norwegischer Mathematiker, hat 1892 als 18jähriger jene Formel für die Zahl Pi gefunden, die dem aktuellen Rekordhalter Prof. Yasumasa Kanada (Univ. Tokyo) 2002 erlaubte, Pi auf 1,241,100,000,000 Dezimalstellen genau zu berechnen :s2:. Størmers Formel lautet
    Pi = 176*arctan(1/57)+ 28*arctan(1/239)-48*arctan(1/682)+96*arctan(1/12943)
    Bei der Auswertung wird üblicherweise die Taylor-Entwicklung der arctan()-Funktion verwendet:
    arctan(x) = x - x^3/3 + x^5/5 - x^7/7 + ...
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3
  • 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 des 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/Randbedingungen 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

    Dieser Beitrag wurde bereits 2 mal editiert, zuletzt von p. specht ()

  • Abt. Datenstruktur DEQue: Nachtrag DEQ.inc Include
    =================================

    Quellcode

    1. Proc DEQHelp
    2. font 2:Cls:Print "\n\n============================================================"
    3. Print " D E Q H e l p"
    4. Print "============================================================"
    5. Print "\n CLEAR q&[] ' natives XProfan ":Print " sz& = SizeOf(q&[]) ' sz=5 bedeutet: 0,1,2,3,4 "
    6. Print " SetSize q&[],sizeof(q&[])+1 ' erweitert Array um 1 Stelle "
    7. Print " MakeIndexDEQ(q&[],n&) ":Print " MakeRevertIndexDEQ(q&[],n&) "
    8. Print " MakeConstDEQ(q&[],n&,wert&) ":Print "\n ShowDEQ q&[],nw& ' Temporary NumWidth "
    9. Print "\n Wert&= PopBtm(q&[]) ' -999999999 = OutOfRange Error "
    10. Print " Wert&= PopTop(q&[]) ' or -999999999 ":Print " PutTop Wert&,q&[] "
    11. Print " PutBtm Wert&,q&[] ":Print "\n Wert&= CopyBtm(q&[]) "
    12. Print " Wert&= CopyTop(q&[]) ' -999999999 ":Print " Wert&= CopyFromDEQ(q&[],n&) "
    13. Print "\n LogVgl&=CmpInDEQ(q&[],n&,m&)' n=Pos1,m=Pos2: <1,'=0',>-1 "
    14. Print "\n NisM&= SwapInDEQ(q&[],n&,m&)' (N=M)=1, -999999999 "
    15. Print " Wert&= DelFromDEQ(q&[],n&) ":Print " DEQHelp ' diesen Text ausgeben":waitinput :endproc
    16. :proc MakeIndexDEQ :parameters q&[],n&:SetSize q&[],n&:q&[]=&index+1:endproc
    17. :proc MakeRevertIndexDEQ :parameters q&[],n&:SetSize q&[],n&:q&[]=n&-&index:endproc
    18. :proc MakeConstDEQ :parameters q&[],n&,wert&:SetSize q&[],n&:q&[]=Wert&:endproc
    19. :Proc PopBtm :parameters q&[]:var sz&=sizeof(q&[]):if sz&>0:var wert&=q&[0]
    20. q&[]=q&[&index+(&index<(sz&-1))]:setsize q&[],sizeof(q&[])-(sz&>0):return wert&
    21. else :return int(-999999999):endif :endproc
    22. :proc PopTop :parameters q&[]:declare wert&:var sz&=sizeof(q&[]):if sz&>0:wert&=q&[sz&-1]
    23. else :wert&=int(-999999999):endif :setsize q&[],sz&-(sz&>0):return wert&:endproc
    24. :Proc PutTop :parameters Wert&,q&[]:var sz&=sizeof(q&[]):inc sz&:setsize q&[],sz&:q&[sz&-1]=Wert&:endproc
    25. :Proc PutBtm :parameters Wert&,q&[]:declare tmp&[]:var sz&=sizeof(q&[]):inc sz&
    26. setsize q&[],sz&:tmp&[]=q&[]:q&[]=tmp&[&index-(&index>0)]:clear tmp&[]:q&[0]=Wert&:endproc
    27. :Proc CopyBtm :parameters q&[]:var sz&=sizeof(q&[]):if sz&>0:return q&[0]
    28. else :return int(-999999999):endif :endproc
    29. :proc CopyTop :parameters q&[]:declare wert&:var sz&=sizeof(q&[])
    30. if sz&>0:wert&=q&[sz&-1]:else :wert&=int(-999999999):endif :return wert&:endproc
    31. :proc ShowDEQ :parameters q&[],nw&:var tmp&=get("numwidth"):set("numwidth",nw&)
    32. var sz&=sizeof(q&[]):print tab(2);"n=";sz&;": ";:whileloop 0,sz&-1:print q&[&Loop],
    33. endwhile :print:set("numwidth",tmp&):endproc
    34. :proc CopyFromDEQ :parameters q&[],n&:declare wert&:var sz&=sizeof(q&[])
    35. if (sz&>0) and (n&>=0) and (n&<sz&):wert&=q&[n&]:else :wert&=int(-999999999):endif :return wert&:endproc
    36. :Proc CmpInDEQ :parameters q&[],n&,m&
    37. var sz&=sizeof(q&[]):if (sz&>0) and (n&>=0) and (n&<sz&) and (m&>=0) and (m&<sz&)
    38. return int((q&[n&]<q&[m&])-(q&[n&]>q&[m&])):else :return int(-999999999):endif :Endproc
    39. :Proc SwapInDEQ :parameters q&[],n&,m& 'Swaps two elements in Queue
    40. var sz&=sizeof(q&[]):if (sz&>0) and (n&>=0) and (n&<sz&) and (m&>=0) and (m&<sz&):var tmp&=q&[n&]
    41. q&[n&]=q&[m&]:q&[m&]=tmp&:return int(n&<>m&):else :return int(-999999999):endif :endproc
    42. :proc DelFromDEQ :parameters q&[],n&:var sz&=sizeof(q&[]) :declare wert&
    43. if (sz&>0) and (n&>0) and (n&<sz&):wert&=q&[n&]:q&[]=q&[&index+(&index>=n&)*(&index<(sz&-1))]
    44. dec sz&:setsize q&[],sz&:else :wert&=int(-999999999):endif :return wert&
    45. endproc
    Alles anzeigen
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3

    Dieser Beitrag wurde bereits 1 mal editiert, zuletzt von Frabbing () aus folgendem Grund: Codetags gesetzt, anstatt Quotetags.

  • Abt. Diverse Links
    ===========
    Html/pdf: Albertson/Hutchinson: Discrete mathematical Algorithms (downloadable version), engl.
    PDF: Prof. Ulrich Hoffmann, Uni Lüneburg: Mathematik für Wirtschaftsinformatiker, deutsch, 235 Seiten
    PDF-Slides: Richard Mayr, Univ. Edinburgh: Discrete Mathematics, Chapter 3: Algorithms, engl., 28 Slides
    PDF: Conrado Martínez, Univ. Politècnica de Catalunya, Spain: Applications of Discrete Mathematics to the Analysis of Algorithms, engl., 40 pages
    PDF-Slides: Prof. Steven Evans, Discrete Mathematics: Number Theory and Cryptography, Chapter Modular Arithmetics, engl.,
    Html mit PDF-Links: Edward A. Bender & S. Gill Williamson: Foundations of Combinatorics with Applications
    Metalink: Free Online Mathematics Books
    ACM TOMS Algorithm 419: Jenkins-Traub-Algorithm Polynomial Root Solver (Free Windows Version with setup.exe)
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3

    Dieser Beitrag wurde bereits 3 mal editiert, zuletzt von p. specht ()

  • P.S. zu oben:
    ========
    PDF-Slides: Dominic Schuhmacher et.al., Univ. Zürich: Die Irrfahrt - eine Entdeckungsreise in die Welt des Zufalls, deutsch, 75 Folien (U.a.: Wann zerfließt ein Wassertropfen?)
    PDF: Stefanie Endres (Hausarbeit, Univ. Würzburg): Irrfahrten: Alle Wege führen nach Rom, deutsch, 114 Seiten, mit zahlreichen Beispielen
    PDF-Slides: Beatrice Paoli, FWS-Referat: Methoden zur Analyse sozialer Netzwerke, deutsch, 42 Folien, (Impfstrategien, Mutierende Zufallsnetzwerke)
    Wikipedia: Petrinetze
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3
  • Abt. Test-StringArray$
    ==============
    Manchmal braucht man mal schnell 15000 Teststrings zufälliger Länge...

    Quellcode

    1. CLS:randomize:font 2
    2. Declare a$[14999]
    3. a$[]=mkstr$(chr$(32+rnd(128)),rnd(55))
    4. a$[]=a$[]+printme():waitinput
    5. :proc printme :print a$[&index]:endproc

    Das wievielte Array (Laufnummer in Reihenfolge wie Declared) ist d![] ?:

    Quellcode

    1. cls:declare a%[],b$,c&[],d![],e$[10],f$[]
    2. var ArrayNr&=d![]:print " d![] ist als ";ArrayNr&;". Array deklariert!":waitinput
    Gruss

    P.S.: Bin auf der Suche nach den Bereichsgrenzen für die Variablendefinitionen, Arrays-Deskriptoren etc. Weiss jemand genaueres? (Es wird wohl heute nicht mehr so funktionieren wie beim Commodore PET oder C=64 von dunnemals, oder? Da standen diese Grenzen sogar im Quick guide/User manual!)
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3

    Dieser Beitrag wurde bereits 4 mal editiert, zuletzt von p. specht ()

  • Und weil Bugs halt kommen und gehen wie sie wollen: Statt a$[]=a$[]+printme():waitinput muss es, um die Inhalte von a$[] zu erhalten, natürlich lauten: a$[]=a$[ &index ]+printme():waitinput
    PrintMe() sollte dann eigentlich auch einen Leerstring zurückliefern, oder man erspart sich die Stringaddition und liefert gleich a$[&index] zurück. Einer der ersten "Bugs" war übrigens keiner, sondern eine Motte - eigentlich müsste es also Moth heissen. "Bug" im Sinne von 'Kleiner Krabbler' (Unberechenbare Kleinigkeit) geht aber angeblich auf Edison und seinen Phonographen zurück.
    Gruss
    Ach so, ja: Abt. Links für Lärm-Liebhaber: Virtual Online-Synthesizer Keyboard. Macht Spaß!
    [IMG:http://upload.wikimedia.org/wikipedia/commons/thumb/8/8a/H96566k.jpg/640px-H96566k.jpg]
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3

    Dieser Beitrag wurde bereits 4 mal editiert, zuletzt von p. specht ()

  • Abt. Fixe String-Arrays: So weit, so breit
    =========================
    Die bisherigen Erkenntnisse betreffend Speicherverwaltung von "mit fester Anzahl deklarierte String-Arrays" (hier: in Profan-11) sind weiterhin dürftig:

    1) 'a%=a$[]' liefert für absolut ALLE Arrays eine Platzziffer in Reihenfolge der declare-ten Arrays.
    2) Der Befehl 'a%=Addr(a$[])' stürzt bei fest declare-ten Stringarrays zwar nicht ab, verweist aber auf Verwaltungsspeicher, der mit Profan-Mitteln nicht gültig auslesbar ist (Exeption Error).
    3) Erkennbar war lediglich, dass jedem String eines 'Feste-Anzahl-Arrays' (Speicherort mit 'a%=Addr(a$[&Loop])' ermittelt) ein Long mit der Länge des Strings VORANGEHT und ...
    4) eine einzelne Null als Endbyte dient (wie in der Hilfe beschrieben).
    5) Ein weiteres Long unmittelbar VOR der Längenangabe enthält stets '1': Vermutlich eine Datentyp-Angabe.

    Gruss

    P.S: Bei dynamischen Stringarrays führt schon das genannte 'a%=Addr(a$[])' zum Absturz mit Exeption error!

    Quellcode

    1. WindowTitle "Fixed-StringArray Memory Testfacility"
    2. Cls:font 0:Randomize:set("randseed",1)
    3. Declare a%,a&,b%[],a$[9],flg%
    4. a$[]=left$(\
    5. chr$(56+rnd(32))+chr$(56+rnd(32))+chr$(56+rnd(32))+chr$(56+rnd(32))\
    6. +chr$(56+rnd(32))+chr$(56+rnd(32))+chr$(56+rnd(32))+chr$(56+rnd(32))\
    7. +chr$(56+rnd(32))+chr$(56+rnd(32))+chr$(56+rnd(32))+chr$(56+rnd(32))\
    8. +chr$(56+rnd(32))+chr$(56+rnd(32)),3+&index)
    9. proc prt
    10. case flg%:print tab(5);&index;".",a$[&index]
    11. inc flg%:return a$[&index]
    12. endproc
    13. flg%=0:a$[]=prt()
    14. flg%=0:a$[]=prt()
    15. a%=a$[]:font 2:print " Das an",a%;". Stelle deklarierte Array!":font 0
    16. whileloop 0,9
    17. a%=addr(a$[&Loop]):print " ";a%,"len:";long(a%,-4),string$(a%,0),"Typ:";long(a%,-8),
    18. print "Endbytes:",byte(addr(a$[&Loop])+long(a%,-4),0);",";byte(addr(a$[&Loop])+long(a%,-4)+1,0);
    19. print ",";byte(addr(a$[&Loop])+long(a%,-4)+2,0);",";byte(addr(a$[&Loop])+long(a%,-4)+3,0)
    20. endwhile
    21. waitinput
    22. end
    Alles anzeigen
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3
  • Abt. Sehenswertes
    ---------------------
    Youtube: 15 Sort-Algorithmen, multimedial in 6 Minuten dargestellt
    Youtube: Quicksort mit Lego in deutsch erklärt
    Html (engl.): The Basics of Game Physics Engines
    Gruss

    P.S.: Es scheint große Mode zu werden, die Soundtracks von Wayne Lyttles 'Animusic' (Computeranimierte Musik) in normalen Orchestern oder Bands nachzuspielen: Beispiel 1 (Amateurband), Beispiel 2 (Orchestral)
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3
  • Abt. Velocity Verlet Integrator
    ==================
    Das große Dilemma bei der numerischen Lösung von Feldgleichungen ist folgendes: Man muss entweder zuerst 2 Orte und die verstrichene Zeit kennen, um eine Durchschnittsgeschwindigkeit berechnen zu können (von der dann aber z.B. die Induktion und auf den Weg zurückwirkende magnetische Kräfte abhängen). Oder man muss bei Kenntnis des Ausgangsortes, der Richtung und Geschwindigkeit den nächsten Bahnpunkt berechnen, die wirkenden nichtlinearen Reibungskräfte sind aber z.B. vom Weg ZWISCHEN Ausgangs- und nächstem Bahnpunkt abhängig ... den man aber noch nicht kennt! Problem erkannt? Fehler-Minimierung allein hilft da wenig, es geht nämlich noch um weitere Eigenschaften, die so ein Simulationsalgorithmus erfüllen muss:

    Das Geheimnis vieler Physik-Simulationen ist Velocity Verlet Integration, ein schnelles Integrationsverfahren für Newton'sche Bewegungsgleichungen, das sogar das hochgenaue 'Runge-Kutta-4./5.Ordnung (RKF45)'-Verfahren übertrifft, - allerdings NICHT an Genauigkeit, sondern an Energie-Konstanz! Runge-Kutta neigt nämlich im Laufe der Zeit dazu, durch systematische Fehler und Rundungsfehler immer mehr Energie in das simulierte System (etwa ein Pendel) zu pumpen.

    Das kann beim vorgestellten VELOCITY VERLET-Verfahren (sehr ähnlich dem 'Leapfrog'-Verfahren, zu deutsch etwa 'Bockspringen') nicht passieren, da es ein ZEITUMKEHRBARES, SYMPLEKTISCHES Verfahren ist, dessen Fehler sich bis auf einen kleinen Anteil gegenseitig weitestgehend aufheben. Zeitumkehrbar heißt: Durch Vorzeichenänderung im Zeitschritt gelangt man wieder an die Ausgangsposition. Symplektisch bedeutet, daß physikalische Erhaltungsgrößen wie Impuls, Energie, Entropie (weitesgehend) erhalten bleiben.

    Gruss
    P.S.: Anbei ein bewusst äußerst einfach gehaltenes Verfahrensbeispiel (In Physik-Engines geht's dann gleich in allen 3 Freiheitsgraden der Translation und 3 der Rotation zur Sache)

    Quellcode

    1. WindowTitle "VELOCITY VERLET DEMO am Beispiel Lineares Feder-Masse-Dämpfungssystem"
    2. '(CL) CopyLeft 2014-08 by P.Specht, Wien
    3. 'Q: http://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet
    4. WindowStyle 24:Window 0,0-%maxx,%maxy
    5. var xh&=width(%hwnd)\2:var yh&=height(%hwnd)\2
    6. Declare m!,dt!,x1!,x2!,v1!,v15!,v2!,a1!,a2!
    7. proc A :parameters x!,v!
    8. var l!=200 ' Federlänge m/1000
    9. var fk!=25 ' Federkonstante N/m 'Kraft pro Auslenkung
    10. var rk!=0.2' Reibungskonstante N*s/m
    11. return (l!-x!)*fk!/m! - v!*rk!
    12. 'Beschleunigung = Kraft / Masse - Reibungsdämpfung
    13. endproc
    14. proc DRAW
    15. usepen 0,(xh&-x2!)/30 ,255
    16. cls
    17. line xh&-500,yh& - xh&+x2!-200,yh&
    18. usepen 0,37,rgb(0,0,255)
    19. line xh&+x2!-200,yh& - xh&+x2!-199,yh&
    20. endproc
    21. 'Startwerte:
    22. dt!=0.03 's
    23. m!=1 'kg
    24. 'Einfachstes Zeichenmodell (mit Cls, ohne MemCopy)
    25. 'Hauptschleife
    26. Repeat :waitinput 33
    27. 'VELOCITY VERLET (aka Leapfrog)
    28. v15!= v1!+a1!/2*dt!
    29. x2! = x1!+v15!*dt!
    30. a2! = A(x2!,v15!)
    31. v2! = v15!+dt!*a2!/2
    32. DRAW
    33. a1!=a2!:v1!=v2!:x1!=x2! ' Umgreifen für nächsten Schritt
    34. until %key=27
    35. END
    Alles anzeigen
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3
  • Abt. AKIMA-Interpolation mit Kubischen Splines
    ------------------------------------------------------
    Ein von US-Volkswirt Hiroshi Akima 1978 vorgestelltes Verfahren nutzt gegenüber der üblichen "linearen Interpolation" zwischen Stützstellen Polynome dritten Grades. In deren Berechnung aus den (auch unregelmäßig verteilten!!) Stützstellen (Messwertetabelle y=F(x)) gehen auch die jeweils benachbarten Stützstellen ein sowie zusätzlich ein gewogener Anstieg, und zwar der Mittelwert aus den Anstiegen zu den benachbarten Stützpunktwerten. Nicht sehr klar? Na, Hauptsache es funktioniert!
    Gruss
    P.S.: In einem revidierten Verfahren hat Akima das 1992 auf Gitterpunkte erweitert, sodaß durch einige diskrete Werte beliebig geformte glatte Flächen definiert werden können, was z.B. im KFZ-Bau verwendet wird.

    Brainfuck-Quellcode

    1. WindowTitle mkstr$(" ",26)+upper$("AKIMA-Interpolation mittels Kubischer Splines [H. Akima, 1978]")
    2. ' Portiert nach XProfan11 (CL) CopyLeft 2014-08 by P. Specht, Vienna (Austria).
    3. '{ Keine wie auch immer geartete Garantie. Nutzung auf alleinige Gefahr des Anwenders!
    4. '
    5. ' Übersetzungsbasis: Turbo Pascal Handbuch 1974; Für Detailklärungen wurde verwendet:
    6. ' F90/95: https://www.tacc.utexas.edu/documents/13601/162125/fortran_class.pdf
    7. ' F95: http://www.mrao.cam.ac.uk/~rachael/compphys/SelfStudyF95.pdf
    8. '
    9. ' Pascal-DemoSource: http://jean-pierre.moreau.pagesperso-orange.fr/Pascal/akima_pas.txt
    10. ' Orig. 2D-Algorithmus: Hiroshi Akima,U.S. Department of Commerce, NTIA/ITS, F90-Version of 1995/08,
    11. ' veröffentlicht in "Hiroshi Akima, 'A Method of Bivariate Interpolation and Smooth Surface Fitting
    12. ' for Irregularly Distributed Data Points", ACM Transactions on Math. Software Vol.4 No.2, June 1978
    13. ' pp. 148-159. Copyright 1978, Association for Computing Machinery, Inc.
    14. ' Fortran-90: TOMS 22.yr,No3 (Sept. 1996) p.357ff, http://calgo.acm.org/760.gz (Abfrg. 2014-08)
    15. ' Test Data: http://cran.r-project.org/web/packages/akima/akima.pdf
    16. '
    17. '********************************************************
    18. '* Program to demonstrate the Akima spline fitting *
    19. '* of Function SIN(X) in double precision *
    20. '* ---------------------------------------------------- *
    21. '* Reference: BASIC Scientific Subroutines, Vol. II *
    22. '* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. *
    23. '* *
    24. '* Pascal Version by J.-P. Moreau, Paris. *
    25. '* (www.jpmoreau.fr) *
    26. '* ---------------------------------------------------- *
    27. '* SAMPLE RUN: *
    28. '* *
    29. '* Akima spline fitting of SIN(X): *
    30. '* *
    31. '* X SIN(X) HANDBOOK AKIMA INTERPOLATION ERROR *
    32. '* ---------------------------------------------------- *
    33. '* 0.00 0.0000000 0.0000000 0.0000000 *
    34. '* 0.05 0.0499792 0.0500402 -0.0000610 *
    35. '* 0.10 0.0998334 0.0998435 -0.0000101 *
    36. '* 0.15 0.1494381 0.1494310 0.0000072 *
    37. '* 0.20 0.1986693 0.1986459 0.0000235 *
    38. '* 0.25 0.2474040 0.2474157 -0.0000118 *
    39. '* 0.30 0.2955202 0.2955218 -0.0000016 *
    40. '* 0.35 0.3428978 0.3428916 0.0000062 *
    41. '* 0.40 0.3894183 0.3894265 -0.0000081 *
    42. '* 0.45 0.4349655 0.4349655 0.0000000 *
    43. '* 0.50 0.4794255 0.4794204 0.0000051 *
    44. '* 0.55 0.5226872 0.5226894 -0.0000021 *
    45. '* 0.60 0.5646425 0.5646493 -0.0000068 *
    46. '* 0.65 0.6051864 0.6051821 0.0000043 *
    47. '* 0.70 0.6442177 0.6442141 0.0000035 *
    48. '* 0.75 0.6816388 0.6816405 -0.0000017 *
    49. '* ---------------------------------------------------- *
    50. '* *
    51. '*******************************************************}
    52. '}
    53. PROC Interpol_Akima
    54. '********************************************************
    55. '* Akima Spline Fitting Subroutine *
    56. '* -----------------------------------------------------*
    57. '* The input table is X(i), Y(i), where Y(i) is the *
    58. '* dependant variable. The interpolation point is xx, *
    59. '* which is assumed to be in the interval of the table *
    60. '* with at least one table value to the left, and three *
    61. '* to the right. The interpolated returned value is yy. *
    62. '* n is returned as an error check (n=0 implies error). *
    63. '* It is also assumed that the X(i) are in ascending *
    64. '* order! *
    65. '********************************************************
    66. parameters i&
    67. n&=1
    68. ' {special case xx=0}
    69. if xx!=0: yy!=0:return:endif
    70. ' {Check to see if interpolation point is correct}
    71. if (xx!<X![1]) or (xx!>=X![iv&-3]): n&=0:return:endif
    72. X![0]=2*X![1]-X![2]
    73. ' {Calculate Akima coefficients, a and b}
    74. whileloop iv&-1:i&=&Loop
    75. ' {Shift i to i+2}
    76. XM![i&+2]=(Y![i&+1]-Y![i&])/(X![i&+1]-X![i&])
    77. endwhile
    78. XM![iv&+2]=2*XM![iv&+1]-XM![iv&]
    79. XM![iv&+3]=2*XM![iv&+2]-XM![iv&+1]
    80. XM![2]=2*XM![3]-XM![4]
    81. XM![1]=2*XM![2]-XM![3]
    82. whileloop iv&:i&=&Loop
    83. a!=ABS(XM![i&+3]-XM![i&+2])
    84. b!=ABS(XM![i&+1]-XM![i&])
    85. if (a!+b!)<>0
    86. Z![i&]=(a!*XM![i&+1]+b!*XM![i&+2])/(a!+b!)
    87. else
    88. Z![i&]=(XM![i&+2]+XM![i&+1])/2
    89. endif
    90. endwhile
    91. ' {Find relevant table interval}
    92. i&=0
    93. repeat
    94. inc i&
    95. until xx!<=X![i&]
    96. dec i&
    97. ' {Begin interpolation}
    98. b!=X![i&+1]-X![i&]
    99. a!=xx!-X![i&]
    100. yy!=Y![i&]+Z![i&]*a!+(3*XM![i&+2]-2*Z![i&]-Z![i&+1])*a!*a!/b!
    101. yy!=yy!+(Z![i&]+Z![i&+1]-2*XM![i&+2])*a!*a!*a!/(b!*b!)
    102. EndProc
    103. '{ Main Program AKIMA-SINTEST
    104. Windowstyle 24:font 0:set("decimals",8)
    105. CLS 'Window 0,0-%maxx,%maxy
    106. var fmt$=" 0.0000000;-0.0000000"
    107. var SIZE& = 14
    108. declare X![SIZE&+1],Y![SIZE&+1],XM![SIZE&+4],Z![SIZE&+1]
    109. declare a!,b!,xx!,yy!, i&,iv&,n&
    110. iv&=14 '{Number of points in table}
    111. CLS
    112. Font 2:Print "\n";tab(10);" Hiroshi Akima's Interpolation mit kubischen Splines"
    113. Font 0:print tab(10);" "+mkstr$("-",53)
    114. print tab(10);" Basis sind 14 unregelmäßig verteilte Stützwerte einer "
    115. print tab(10);" Eigenbau-Sinusfunktion. Errechnet wird der absolute "
    116. Print tab(10);" Fehler (als Differenz zur internen Sin(x)-Funktion). "
    117. print tab(10);" Im Programmtext sind die Vergleichswerte zu finden! \n"
    118. ' -----------------------------------------------------------------
    119. ' Input y=Sine(x) table
    120. ' -----------------------------------------------------------------
    121. ' Sine table values from Handbook of mathematical functions
    122. ' by M. Abramowitz and I.A. Stegun, NBS, june 1964
    123. ' -----------------------------------------------------------------
    124. X![1]=0.000:Y![1]=0.00000000 : X![2]=0.125:Y![2]=0.12467473
    125. X![3]=0.217:Y![3]=0.21530095 : X![4]=0.299:Y![4]=0.29456472
    126. X![5]=0.376:Y![5]=0.36720285 : X![6]=0.450:Y![6]=0.43496553
    127. X![7]=0.520:Y![7]=0.49688014 : X![8]=0.589:Y![8]=0.55552980
    128. X![9]=0.656:Y![9]=0.60995199 : X![10]=0.721:Y![10]=0.66013615
    129. X![11]=0.7853981634 : Y![11]=0.7071067812
    130. X![12]=0.849:Y![12]=0.75062005: X![13]=0.911:Y![13]=0.79011709
    131. X![14]=0.972:Y![14]=0.82601466
    132. ' -----------------------------------------------------------------
    133. print "\n";mkstr$("-",77)
    134. font 2
    135. print tab(3+3);"X";tab(%pos+15);"SIN(X)";\
    136. tab(%pos+3);"AKIMA-Stützwertinterpolation";tab(%pos+2);"ABS.FEHLER"
    137. font 0
    138. print mkstr$("=",77)
    139. ' {main loop}
    140. xx!=0
    141. whileloop 16:i&=&Loop
    142. Interpol_Akima(i&)
    143. print tab(4);xx!;tab(%pos+6);SIN(xx!);tab(%pos+9);yy!;tab(%pos+11);format$(fmt$,SIN(xx!) - yy!)
    144. casenot i& mod 4:print
    145. xx! = xx! + 0.05
    146. endwhile
    147. ' {print footer}
    148. print mkstr$("-",77)
    149. beep
    150. waitinput
    151. END
    152. '} of file akima.prf
    Alles anzeigen
    Win7-64HomPremSP1,XProfan11.2a,XPIA,JWasm,xpse,IntelCoreQuad2.5GHz/4GB/je1TB HD intern:esataBay:USB3