import matplotlib.pyplot as plt

def model():
    """Solve u'' = -1, u(0)=0, u'(1)=0."""
    import sympy as sym
    x, c_0, c_1, = sym.symbols('x c_0 c_1')
    u_x = sym.integrate(1, (x, 0, x)) + c_0
    u = sym.integrate(u_x, (x, 0, x)) + c_1
    r = sym.solve([u.subs(x,0) - 0,
                   sym.diff(u,x).subs(x, 1) - 0],
                  [c_0, c_1])
    u = u.subs(c_0, r[c_0]).subs(c_1, r[c_1])
    u = sym.simplify(sym.expand(u))
    return u

def midpoint_rule(f, M=100000):
    """Integrate f(x) over [0,1] using M intervals."""
    from numpy import sum, linspace
    dx = 1.0/M                       # interval length
    x = linspace(dx/2, 1-dx/2, M)    # integration points
    return dx*sum(f(x))

def check_integral_b():
    from numpy import pi, sin
    for i in range(12):
        exact = 2/(pi*(2*i+1))
        numerical = midpoint_rule(
            f=lambda x: sin((2*i+1)*pi*x/2))
        print(i, abs(exact - numerical))

def sine_sum(x, N):
    s = 0
    from numpy import pi, sin, zeros
    u = [] # u[k] is the sum i=0,...,k
    k = 0
    for i in range(N+1):
        s += - 16.0/((2*i+1)**3*pi**3)*sin((2*i+1)*pi*x/2)
        u.append(s.copy())  # important with copy!
    return u

def plot_sine_sum():
    from numpy import linspace
    x = linspace(0, 1, 501)  # coordinates for plot
    u = sine_sum(x, N=10)
    u_e = 0.5*x*(x-2)
    N_values = 0, 1, 10
    for k in N_values:
        plt.plot(x, u[k])
    plt.plot(x, u_e)
    plt.legend(['N=%d' % k for k in N_values] + ['exact'],
               loc='upper right')
    plt.xlabel('$x$');  plt.ylabel('$u$')
    plt.savefig('tmpc.png'); plt.savefig('tmpc.pdf')

def check_integral_d():
    from numpy import pi, sin
    for i in range(24):
        if i % 2 == 0:
            exact = 2/(pi*(i+1))
        elif (i-1) % 4 == 0:
            exact = 2*2/(pi*(i+1))
        else:
            exact = 0
        numerical = midpoint_rule(
            f=lambda x: sin((i+1)*pi*x/2))
        print(i, abs(exact - numerical))

def check_integral_d_sympy_answer():
    from numpy import pi, sin
    for i in range(12):
        exact = 2/(pi*(i+1))
        numerical = midpoint_rule(
            f=lambda x: sin((i+1)*pi*x/2))
        print(i, abs(exact - numerical))

def sine_sum_d(x, N):
    s = 0
    from numpy import pi, sin, zeros
    u = []  # u[k] is the sum i=0,...,k
    k = 0
    for i in range(N+1):
        if i % 2 == 0:       # even i
            s +=   - 16.0/((i+1)**3*pi**3)*sin((i+1)*pi*x/2)
        elif (i-1) % 4 == 0:   # 1, 5, 9, 13, 17
            s += - 2*16.0/((i+1)**3*pi**3)*sin((i+1)*pi*x/2)
        else:
            s += 0
        u.append(s.copy())
    return u

def plot_sine_sum_d():
    from numpy import linspace
    x = linspace(0, 1, 501)  # coordinates for plot
    u = sine_sum_d(x, N=20)
    u_e = 0.5*x*(x-2)
    N_values = 0, 1, 2, 3, 20
    for k in N_values:
        plt.plot(x, u[k])
    plt.plot(x, u_e)
    plt.legend(['N=%d' % k for k in N_values] + ['exact'],
               loc='upper right')
    plt.xlabel('$x$');  plt.ylabel('$u$')
    #plt.axis([0.9, 1, -0.52, -0.49])
    plt.savefig('tmpd.png'); plt.savefig('tmpd.pdf')


if __name__ == '__main__':
    import sys
    print(model())
    print('sine 2*i+1 integral:')
    check_integral_b()
    print('sine i+1 integral, sympy answer:')
    check_integral_d_sympy_answer()
    print('sine i+1 integral:')
    check_integral_d()
    #sys.exit(0)
    plot_sine_sum()
    plt.figure()
    plot_sine_sum_d()
    plt.show()

[((132, 156), 'sympy.symbols', 'sym.symbols', (['"""x c_0 c_1"""'], {}), "('x c_0 c_1')\n", (143, 156), True, 'import sympy as sym\n'), ((643, 674), 'numpy.linspace', 'linspace', (['(dx / 2)', '(1 - dx / 2)', 'M'], {}), '(dx / 2, 1 - dx / 2, M)\n', (651, 674), False, 'from numpy import linspace\n'), ((1274, 1293), 'numpy.linspace', 'linspace', (['(0)', '(1)', '(501)'], {}), '(0, 1, 501)\n', (1282, 1293), False, 'from numpy import linspace\n'), ((1443, 1459), 'matplotlib.pyplot.plot', 'plt.plot', (['x', 'u_e'], {}), '(x, u_e)\n', (1451, 1459), True, 'import matplotlib.pyplot as plt\n'), ((1464, 1539), 'matplotlib.pyplot.legend', 'plt.legend', (["([('N=%d' % k) for k in N_values] + ['exact'])"], {'loc': '"""upper right"""'}), "([('N=%d' % k) for k in N_values] + ['exact'], loc='upper right')\n", (1474, 1539), True, 'import matplotlib.pyplot as plt\n'), ((1557, 1574), 'matplotlib.pyplot.xlabel', 'plt.xlabel', (['"""$x$"""'], {}), "('$x$')\n", (1567, 1574), True, 'import matplotlib.pyplot as plt\n'), ((1577, 1594), 'matplotlib.pyplot.ylabel', 'plt.ylabel', (['"""$u$"""'], {}), "('$u$')\n", (1587, 1594), True, 'import matplotlib.pyplot as plt\n'), ((1599, 1622), 'matplotlib.pyplot.savefig', 'plt.savefig', (['"""tmpc.png"""'], {}), "('tmpc.png')\n", (1610, 1622), True, 'import matplotlib.pyplot as plt\n'), ((1624, 1647), 'matplotlib.pyplot.savefig', 'plt.savefig', (['"""tmpc.pdf"""'], {}), "('tmpc.pdf')\n", (1635, 1647), True, 'import matplotlib.pyplot as plt\n'), ((2732, 2751), 'numpy.linspace', 'linspace', (['(0)', '(1)', '(501)'], {}), '(0, 1, 501)\n', (2740, 2751), False, 'from numpy import linspace\n'), ((2909, 2925), 'matplotlib.pyplot.plot', 'plt.plot', (['x', 'u_e'], {}), '(x, u_e)\n', (2917, 2925), True, 'import matplotlib.pyplot as plt\n'), ((2930, 3005), 'matplotlib.pyplot.legend', 'plt.legend', (["([('N=%d' % k) for k in N_values] + ['exact'])"], {'loc': '"""upper right"""'}), "([('N=%d' % k) for k in N_values] + ['exact'], loc='upper right')\n", (2940, 3005), True, 'import matplotlib.pyplot as plt\n'), ((3023, 3040), 'matplotlib.pyplot.xlabel', 'plt.xlabel', (['"""$x$"""'], {}), "('$x$')\n", (3033, 3040), True, 'import matplotlib.pyplot as plt\n'), ((3043, 3060), 'matplotlib.pyplot.ylabel', 'plt.ylabel', (['"""$u$"""'], {}), "('$u$')\n", (3053, 3060), True, 'import matplotlib.pyplot as plt\n'), ((3103, 3126), 'matplotlib.pyplot.savefig', 'plt.savefig', (['"""tmpd.png"""'], {}), "('tmpd.png')\n", (3114, 3126), True, 'import matplotlib.pyplot as plt\n'), ((3128, 3151), 'matplotlib.pyplot.savefig', 'plt.savefig', (['"""tmpd.pdf"""'], {}), "('tmpd.pdf')\n", (3139, 3151), True, 'import matplotlib.pyplot as plt\n'), ((3450, 3462), 'matplotlib.pyplot.figure', 'plt.figure', ([], {}), '()\n', (3460, 3462), True, 'import matplotlib.pyplot as plt\n'), ((3489, 3499), 'matplotlib.pyplot.show', 'plt.show', ([], {}), '()\n', (3497, 3499), True, 'import matplotlib.pyplot as plt\n'), ((167, 194), 'sympy.integrate', 'sym.integrate', (['(1)', '(x, 0, x)'], {}), '(1, (x, 0, x))\n', (180, 194), True, 'import sympy as sym\n'), ((209, 238), 'sympy.integrate', 'sym.integrate', (['u_x', '(x, 0, x)'], {}), '(u_x, (x, 0, x))\n', (222, 238), True, 'import sympy as sym\n'), ((428, 441), 'sympy.expand', 'sym.expand', (['u'], {}), '(u)\n', (438, 441), True, 'import sympy as sym\n'), ((1421, 1438), 'matplotlib.pyplot.plot', 'plt.plot', (['x', 'u[k]'], {}), '(x, u[k])\n', (1429, 1438), True, 'import matplotlib.pyplot as plt\n'), ((2887, 2904), 'matplotlib.pyplot.plot', 'plt.plot', (['x', 'u[k]'], {}), '(x, u[k])\n', (2895, 2904), True, 'import matplotlib.pyplot as plt\n'), ((1129, 1158), 'numpy.sin', 'sin', (['((2 * i + 1) * pi * x / 2)'], {}), '((2 * i + 1) * pi * x / 2)\n', (1132, 1158), False, 'from numpy import pi, sin, zeros\n'), ((2468, 2493), 'numpy.sin', 'sin', (['((i + 1) * pi * x / 2)'], {}), '((i + 1) * pi * x / 2)\n', (2471, 2493), False, 'from numpy import pi, sin, zeros\n'), ((886, 915), 'numpy.sin', 'sin', (['((2 * i + 1) * pi * x / 2)'], {}), '((2 * i + 1) * pi * x / 2)\n', (889, 915), False, 'from numpy import pi, sin, zeros\n'), ((1942, 1967), 'numpy.sin', 'sin', (['((i + 1) * pi * x / 2)'], {}), '((i + 1) * pi * x / 2)\n', (1945, 1967), False, 'from numpy import pi, sin, zeros\n'), ((2182, 2207), 'numpy.sin', 'sin', (['((i + 1) * pi * x / 2)'], {}), '((i + 1) * pi * x / 2)\n', (2185, 2207), False, 'from numpy import pi, sin, zeros\n'), ((2578, 2603), 'numpy.sin', 'sin', (['((i + 1) * pi * x / 2)'], {}), '((i + 1) * pi * x / 2)\n', (2581, 2603), False, 'from numpy import pi, sin, zeros\n'), ((300, 314), 'sympy.diff', 'sym.diff', (['u', 'x'], {}), '(u, x)\n', (308, 314), True, 'import sympy as sym\n')]