finite_differences.py 5.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135
  1. """
  2. Illustrate forward, backward and centered finite differences
  3. in four figures.
  4. """
  5. from pysketcher import *
  6. #test_test()
  7. xaxis = 2
  8. drawing_tool.set_coordinate_system(0, 7, 1, 6, axis=False)
  9. f = SketchyFunc1('$u(t)$')
  10. x = 3 # center point where we want the derivative
  11. xb = 2 # x point used for backward difference
  12. xf = 4 # x point used for forward difference
  13. p = (x, f(x)) # center point
  14. pf = (xf, f(xf)) # forward point
  15. pb = (xb, f(xb)) # backward point
  16. r = 0.1 # radius of circles placed at key points
  17. c = Circle(p, r).set_linecolor('blue')
  18. cf = Circle(pf, r).set_linecolor('red')
  19. cb = Circle(pb, r).set_linecolor('green')
  20. # Points in the mesh
  21. p0 = point(x, xaxis) # center point
  22. pf0 = point(xf, xaxis) # forward point
  23. pb0 = point(xb, xaxis) # backward point
  24. tick = 0.05
  25. # 1D mesh with three points
  26. mesh = Composition({
  27. 'tnm1': Text('$t_{n-1}$', pb0 - point(0, 0.3)),
  28. 'tn': Text('$t_{n}$', p0 - point(0, 0.3)),
  29. 'tnp1': Text('$t_{n+1}$', pf0 - point(0, 0.3)),
  30. 'axis': Composition({
  31. 'hline': Line(pf0-point(3,0), pb0+point(3,0)).\
  32. set_linecolor('black').set_linewidth(1),
  33. 'tick_m1': Line(pf0+point(0,tick), pf0-point(0,tick)).\
  34. set_linecolor('black').set_linewidth(1),
  35. 'tick_n': Line(p0+point(0,tick), p0-point(0,tick)).\
  36. set_linecolor('black').set_linewidth(1),
  37. 'tick_p1': Line(pb0+point(0,tick), pb0-point(0,tick)).\
  38. set_linecolor('black').set_linewidth(1)}),
  39. })
  40. # 1D mesh with three points for Crank-Nicolson
  41. mesh_cn = Composition({
  42. 'tnm1': Text('$t_{n}$', pb0 - point(0, 0.3)),
  43. 'tn': Text(r'$t_{n+\frac{1}{2}}$', p0 - point(0, 0.3)),
  44. 'tnp1': Text('$t_{n+1}$', pf0 - point(0, 0.3)),
  45. 'axis': Composition({
  46. 'hline': Line(pf0-point(3,0), pb0+point(3,0)).\
  47. set_linecolor('black').set_linewidth(1),
  48. 'tick_m1': Line(pf0+point(0,tick), pf0-point(0,tick)).\
  49. set_linecolor('black').set_linewidth(1),
  50. 'tick_n': Line(p0+point(0,tick), p0-point(0,tick)).\
  51. set_linecolor('black').set_linewidth(1),
  52. 'tick_p1': Line(pb0+point(0,tick), pb0-point(0,tick)).\
  53. set_linecolor('black').set_linewidth(1)}),
  54. })
  55. # Vertical dotted lines at each mesh point
  56. vlinec = Line(p, p0).set_linestyle('dotted').\
  57. set_linecolor('blue').set_linewidth(1)
  58. vlinef = Line(pf, pf0).set_linestyle('dotted').\
  59. set_linecolor('red').set_linewidth(1)
  60. vlineb = Line(pb, pb0).set_linestyle('dotted').\
  61. set_linecolor('green').set_linewidth(1)
  62. # Compose vertical lines for each type of difference
  63. forward_lines = Composition({'center': vlinec, 'right': vlinef})
  64. backward_lines = Composition({'center': vlinec, 'left': vlineb})
  65. centered_lines = Composition({'left': vlineb, 'right': vlinef})
  66. centered_lines2 = Composition({'left': vlineb, 'right': vlinef,
  67. 'center': vlinec})
  68. # Tangents illustrating the derivative
  69. domain = [1, 5]
  70. domain2 = [2, 5]
  71. forward_tangent = Line(p, pf).new_interval(x=domain2).\
  72. set_linestyle('dashed').set_linecolor('red')
  73. backward_tangent = Line(pb, p).new_interval(x=domain).\
  74. set_linestyle('dashed').set_linecolor('green')
  75. centered_tangent = Line(pb, pf).new_interval(x=domain).\
  76. set_linestyle('dashed').set_linecolor('blue')
  77. h = 1E-3 # h in finite difference approx used to compute the exact tangent
  78. exact_tangent = Line((x+h, f(x+h)), (x-h, f(x-h))).\
  79. new_interval(x=domain).\
  80. set_linestyle('dotted').set_linecolor('black')
  81. forward = Composition(
  82. dict(tangent=forward_tangent,
  83. point1=c, point2=cf, coor=forward_lines,
  84. name=Text('forward',
  85. forward_tangent.geometric_features()['end'] + \
  86. point(0.1,0), alignment='left')))
  87. backward = Composition(
  88. dict(tangent=backward_tangent,
  89. point1=c, point2=cb, coor=backward_lines,
  90. name=Text('backward',
  91. backward_tangent.geometric_features()['end'] + \
  92. point(0.1,0), alignment='left')))
  93. centered = Composition(
  94. dict(tangent=centered_tangent,
  95. point1=cb, point2=cf, point=c, coor=centered_lines2,
  96. name=Text('centered',
  97. centered_tangent.geometric_features()['end'] + \
  98. point(0.1,0), alignment='left')))
  99. exact = Composition(dict(graph=f, tangent=exact_tangent))
  100. forward = Composition(dict(difference=forward, exact=exact)).\
  101. set_name('forward')
  102. backward = Composition(dict(difference=backward, exact=exact)).\
  103. set_name('backward')
  104. centered = Composition(dict(difference=centered, exact=exact)).\
  105. set_name('centered')
  106. all = Composition(
  107. dict(exact=exact, forward=forward, backward=backward,
  108. centered=centered)).set_name('all')
  109. for fig in forward, backward, centered, all:
  110. drawing_tool.erase()
  111. fig.draw()
  112. mesh.draw()
  113. drawing_tool.display()
  114. drawing_tool.savefig('fd_'+fig.get_name())
  115. # Crank-Nicolson around t_n+1/2
  116. drawing_tool.erase()
  117. centered.draw()
  118. mesh_cn.draw()
  119. drawing_tool.display()
  120. drawing_tool.savefig('fd_centered_CN')
  121. raw_input()