Hans Petter Langtangen 10 年之前
父節點
當前提交
5be557eebb

二進制
examples/FE_comic_strip.pdf


二進制
examples/FE_comic_strip.png


二進制
examples/FE_noncomic_strip.pdf


二進制
examples/FE_noncomic_strip.png


+ 83 - 0
examples/ForwardEuler.py

@@ -0,0 +1,83 @@
+from pysketcher import *
+import numpy as np
+
+xmin = 0
+drawing_tool.set_coordinate_system(xmin=xmin, xmax=4,
+                                   ymin=0, ymax=2.5,
+                                   axis=True, xkcd=True)
+drawing_tool.set_linecolor('blue')
+
+def ForwardEuler(I, a, T, dt):
+    u = [I]
+    t = [0]
+    while t[-1] <= T:
+        u_new = u[-1] - a*dt*u[-1]
+        u.append(u_new)
+        t.append(t[-1] + dt)
+    return np.array(u), np.array(t)
+
+def make_fig(dt=0.5):
+    I = 2
+    a = 0.5
+    T_e = 3
+    T_FE = 1
+    t_fine = np.linspace(0, T_e, 101)
+    u_e = I*np.exp(-a*t_fine)
+
+    u, t = ForwardEuler(I, a, T_FE, dt)
+
+    # y = slope*(x - x0) + y0
+    # u_extrapolated = -a*u[-1]*(t - t[-1]) + u[-1]
+    t_future = t[-1] + 1.5   # let the line be longer than one step
+    line = Line((t[-1], u[-1]), (t_future, -a*u[-1]*(t_future - t[-1]) + u[-1]))
+
+    circles = {
+        i: Circle((t[i], u[i]), 0.05).set_linecolor('red').set_filled_curves('red')
+        for i in range(1, len(u))}
+    # Add next predicted point
+    t_next = t[-1] + dt
+    u_next = -a*u[-1]*(t_next - t[-1]) + u[-1]
+    circles[0] = Circle((t_next, u_next), 0.05).\
+                 set_linecolor('red').set_filled_curves('red')
+    circles = Composition(circles)
+
+    curves = Composition(dict(
+        exact=Curve(t_fine, u_e).set_linestyle('dashed'),
+        numerical=Curve(t, u),
+        extrapolation=line.set_linecolor('red').set_linewidth(3)))
+
+    text_exact = Text_wArrow("exact solution", (2.5, 1), (2.5, I*np.exp(-a*2.5)),
+                             alignment='left')
+
+    text_predict = Text_wArrow("Here we know the slope:\n$u'=-au$!\nLet the solution continue\nalong that slope.",
+                               (1.7, 1.7), (t[-1], u[-1]),
+                               alignment='left')
+    text_next = Text_wArrow("This is the next\npredicted point",
+                            (1, 0.25), (t_next, u_next),
+                            alignment='left')
+
+    fig = Composition(dict(curves=curves,
+                           circles=circles,
+                           exact=text_exact,
+                           predict=text_predict,
+                           next=text_next))
+    return fig
+
+fig = make_fig(dt=0.5)
+fig.draw()
+drawing_tool.display()
+drawing_tool.savefig('tmp1')
+
+drawing_tool.erase()
+text_comment = Text('Just reduce the time step to\nmake more accurate predictions!', (1, 2.25))
+fig = make_fig(dt=0.24)
+fig = Composition(dict(prev_fig=fig, comment=text_comment))
+fig.draw()
+drawing_tool.display()
+drawing_tool.savefig('tmp2')
+
+import os
+os.system('doconce combine_images pdf -2 tmp1 tmp2 FE_strip')
+os.system('doconce combine_images png -2 tmp1 tmp2 FE_strip')
+
+raw_input()

+ 11 - 11
examples/beam2.py

@@ -45,22 +45,22 @@ M1.set_linecolor('black')
 ab_level = point(0, 3*h)
 a_dim = Distance_wText(A - ab_level, B - ab_level, '$a$')
 b_dim = Distance_wText(B - ab_level, C - ab_level, '$b$')
-dims = Compose({'a': a_dim, 'b': b_dim})
-symbols = Compose({'R1': R1, 'R2': R2, 'M1': M1,
-                   'w': load, 'w text': load_text,
-                   'A': Text('$A$', A+point(0.7*h,-0.9*h)),
-                   'B': Text('$B$', support.mid_support-point(h,0)),
-                   'C': Text('$C$', C+point(h/2,-h/2))})
+dims = Composition({'a': a_dim, 'b': b_dim})
+symbols = Composition({'R1': R1, 'R2': R2, 'M1': M1,
+                       'w': load, 'w text': load_text,
+                       'A': Text('$A$', A+point(0.7*h,-0.9*h)),
+                       'B': Text('$B$', support.mid_support-point(h,0)),
+                       'C': Text('$C$', C+point(h/2,-h/2))})
 
 x_axis = Axis(A + point(L+h, H/2), 2*H, '$x$',).set_linecolor('black')
 y_axis = Axis(A + point(0,H/2), 3.5*H, '$y$',
               below=False, rotation_angle=90).set_linecolor('black')
-axes = Compose({'x axis': x_axis, 'y axis': y_axis})
+axes = Composition({'x axis': x_axis, 'y axis': y_axis})
 
-annotations = Compose({'dims': dims, 'symbols': symbols,
-                'axes': axes})
-fig = Compose({'beam': beam, 'support': support,
-               'clamped end': clamped, 'load': load})
+annotations = Composition({'dims': dims, 'symbols': symbols,
+                           'axes': axes})
+fig = Composition({'beam': beam, 'support': support,
+                   'clamped end': clamped, 'load': load})
 
 def deflection(x, a, b, w):
     import numpy as np

+ 83 - 0
examples/integral.py

@@ -0,0 +1,83 @@
+from pysketcher import *
+
+def f(x):
+    return 3*np.exp(-x**4)
+
+xmin = -2
+drawing_tool.set_coordinate_system(xmin=xmin, xmax=4,
+                                   ymin=0, ymax=4,
+                                   axis=True, xkcd=True)
+drawing_tool.set_linecolor('blue')
+
+import numpy as np
+x = np.linspace(xmin, 3, 201)
+y = f(x)
+curve = Curve(x, y)
+
+x2 = np.linspace(xmin, 0.2, 201)
+y2 = f(x2)
+x2 = x2.tolist()
+y2 = y2.tolist()
+x2.append(x2[-1])
+y2.append(0)
+x2.append(xmin)
+y2.append(0)
+#x2 = np.array(x2)
+#y2 = np.array(y2)
+filled = Curve(x2, y2).set_filled_curves(pattern='/')
+
+text1 = Text_wArrow('The integral $\int_{-\infty}^{0.2} 3e^{-x^4}dx$\nis impossible to calculate\nby hand but so easy with\na program!', (1.5, 3.5), (-0.2, 1), alignment='left')
+
+fig = Composition(dict(curve=curve, integral=filled, comment=text1))
+fig.draw()
+drawing_tool.display()
+drawing_tool.savefig('tmp1')
+
+#raw_input()
+
+# Draw piecewise curve for midpoint rule
+def piecewise_curve_for_midpoint_rule(N):
+    dx = (0.2 - xmin)/float(N)
+    x3_double = []
+    y3_double = []
+    y_prev = 0
+    for i in range(N):
+        x = xmin + i*dx
+        x3_double.append(x)
+        y3_double.append(y_prev)
+        x3_double.append(x)
+        y_next = 0.5*(f(x) + f(xmin+(i+1)*dx))
+        y3_double.append(y_next)
+        y_prev = y_next
+    x = xmin + (i+1)*dx
+    x3_double.append(x)
+    y3_double.append(y_prev)
+    x3_double.append(x)
+    y3_double.append(0)
+    # Back to start
+    x3_double.append(xmin)
+    y3_double.append(0)
+    midpoint_curve = Curve(x3_double, y3_double).set_filled_curves(pattern='/').set_linecolor('red')
+    return midpoint_curve
+
+text2 = Text_wArrow('We just draw some rectangles\nto approximate the area\nunder the curve and sum up\nthe rectangular areas!', (1.2, 3.5), (-0.2, 1), alignment='left')
+
+drawing_tool.erase()
+fig = Composition(dict(curve=curve, integral=piecewise_curve_for_midpoint_rule(N=4), comment=text2))
+fig.draw()
+drawing_tool.display()
+drawing_tool.savefig('tmp2')
+
+text3 = Text_wArrow('Just add more rectangles\ngo get the integral\nmore accurate!', (1.5, 3.5), (-0.2, 1), alignment='left')
+
+drawing_tool.erase()
+fig = Composition(dict(curve=curve, integral=piecewise_curve_for_midpoint_rule(N=10), comment=text3))
+fig.draw()
+drawing_tool.display()
+drawing_tool.savefig('tmp3')
+
+import os
+os.system('doconce combine_images pdf -3 tmp1 tmp2 tmp3 integral_comic_strip')
+os.system('doconce combine_images png -3 tmp1 tmp2 tmp3 integral_comic_strip')
+
+raw_input()

二進制
examples/integral_comic_strip.pdf


二進制
examples/integral_comic_strip.png


二進制
examples/integral_noncomic_strip.pdf


二進制
examples/integral_noncomic_strip.png


+ 5 - 1
pysketcher/MatplotlibDraw.py

@@ -57,7 +57,8 @@ class MatplotlibDraw:
         self.ax.set_ylim(minmax['ymin']-y_space/2., minmax['ymax']+y_space/2.)
 
     def set_coordinate_system(self, xmin, xmax, ymin, ymax, axis=False,
-                              instruction_file=None, new_figure=True):
+                              instruction_file=None, new_figure=True,
+                              xkcd=False):
         """
         Define the drawing area [xmin,xmax]x[ymin,ymax].
         axis: None or False means that axes with tickmarks
@@ -76,6 +77,9 @@ class MatplotlibDraw:
                 self.instruction_file.close()  # make new py file for commands
 
         self.mpl = mpl
+        if xkcd:
+            self.mpl.xkcd()
+
         self.xmin, self.xmax, self.ymin, self.ymax = \
              float(xmin), float(xmax), float(ymin), float(ymax)
         self.xrange = self.xmax - self.xmin