ForwardEuler.py 2.6 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283
  1. from pysketcher import *
  2. import numpy as np
  3. xmin = 0
  4. drawing_tool.set_coordinate_system(xmin=xmin, xmax=4,
  5. ymin=0, ymax=2.5,
  6. axis=True, xkcd=True)
  7. drawing_tool.set_linecolor('blue')
  8. def ForwardEuler(I, a, T, dt):
  9. u = [I]
  10. t = [0]
  11. while t[-1] <= T:
  12. u_new = u[-1] - a*dt*u[-1]
  13. u.append(u_new)
  14. t.append(t[-1] + dt)
  15. return np.array(u), np.array(t)
  16. def make_fig(dt=0.5):
  17. I = 2
  18. a = 0.5
  19. T_e = 3
  20. T_FE = 1
  21. t_fine = np.linspace(0, T_e, 101)
  22. u_e = I*np.exp(-a*t_fine)
  23. u, t = ForwardEuler(I, a, T_FE, dt)
  24. # y = slope*(x - x0) + y0
  25. # u_extrapolated = -a*u[-1]*(t - t[-1]) + u[-1]
  26. t_future = t[-1] + 1.5 # let the line be longer than one step
  27. line = Line((t[-1], u[-1]), (t_future, -a*u[-1]*(t_future - t[-1]) + u[-1]))
  28. circles = {
  29. i: Circle((t[i], u[i]), 0.05).set_linecolor('red').set_filled_curves('red')
  30. for i in range(1, len(u))}
  31. # Add next predicted point
  32. t_next = t[-1] + dt
  33. u_next = -a*u[-1]*(t_next - t[-1]) + u[-1]
  34. circles[0] = Circle((t_next, u_next), 0.05).\
  35. set_linecolor('red').set_filled_curves('red')
  36. circles = Composition(circles)
  37. curves = Composition(dict(
  38. exact=Curve(t_fine, u_e).set_linestyle('dashed'),
  39. numerical=Curve(t, u),
  40. extrapolation=line.set_linecolor('red').set_linewidth(3)))
  41. text_exact = Text_wArrow("exact solution", (2.5, 1), (2.5, I*np.exp(-a*2.5)),
  42. alignment='left')
  43. text_predict = Text_wArrow("Here we know the slope:\n$u'=-au$!\nLet the solution continue\nalong that slope.",
  44. (1.7, 1.7), (t[-1], u[-1]),
  45. alignment='left')
  46. text_next = Text_wArrow("This is the next\npredicted point",
  47. (1, 0.25), (t_next, u_next),
  48. alignment='left')
  49. fig = Composition(dict(curves=curves,
  50. circles=circles,
  51. exact=text_exact,
  52. predict=text_predict,
  53. next=text_next))
  54. return fig
  55. fig = make_fig(dt=0.5)
  56. fig.draw()
  57. drawing_tool.display()
  58. drawing_tool.savefig('tmp1')
  59. drawing_tool.erase()
  60. text_comment = Text('Just reduce the time step to\nmake more accurate predictions!', (1, 2.25))
  61. fig = make_fig(dt=0.24)
  62. fig = Composition(dict(prev_fig=fig, comment=text_comment))
  63. fig.draw()
  64. drawing_tool.display()
  65. drawing_tool.savefig('tmp2')
  66. import os
  67. os.system('doconce combine_images pdf -2 tmp1 tmp2 FE_strip')
  68. os.system('doconce combine_images png -2 tmp1 tmp2 FE_strip')
  69. raw_input()