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()

[((6, 19, 6, 43), 'sympy.symbols', 'sym.symbols', ({(6, 31, 6, 42): '"""x c_0 c_1"""'}, {}), "('x c_0 c_1')", True, 'import sympy as sym\n'), ((20, 8, 20, 33), 'numpy.linspace', 'linspace', ({(20, 17, 20, 21): 'dx / 2', (20, 23, 20, 29): '1 - dx / 2', (20, 31, 20, 32): 'M'}, {}), '(dx / 2, 1 - dx / 2, M)', False, 'from numpy import linspace\n'), ((43, 8, 43, 27), 'numpy.linspace', 'linspace', ({(43, 17, 43, 18): '0', (43, 20, 43, 21): '1', (43, 23, 43, 26): '501'}, {}), '(0, 1, 501)', False, 'from numpy import linspace\n'), ((49, 4, 49, 20), 'matplotlib.pyplot.plot', 'plt.plot', ({(49, 13, 49, 14): 'x', (49, 16, 49, 19): 'u_e'}, {}), '(x, u_e)', True, 'import matplotlib.pyplot as plt\n'), ((50, 4, 51, 33), 'matplotlib.pyplot.legend', 'plt.legend', (), '', True, 'import matplotlib.pyplot as plt\n'), ((52, 4, 52, 21), 'matplotlib.pyplot.xlabel', 'plt.xlabel', ({(52, 15, 52, 20): '"""$x$"""'}, {}), "('$x$')", True, 'import matplotlib.pyplot as plt\n'), ((52, 24, 52, 41), 'matplotlib.pyplot.ylabel', 'plt.ylabel', ({(52, 35, 52, 40): '"""$u$"""'}, {}), "('$u$')", True, 'import matplotlib.pyplot as plt\n'), ((53, 4, 53, 27), 'matplotlib.pyplot.savefig', 'plt.savefig', ({(53, 16, 53, 26): '"""tmpc.png"""'}, {}), "('tmpc.png')", True, 'import matplotlib.pyplot as plt\n'), ((53, 29, 53, 52), 'matplotlib.pyplot.savefig', 'plt.savefig', ({(53, 41, 53, 51): '"""tmpc.pdf"""'}, {}), "('tmpc.pdf')", True, 'import matplotlib.pyplot as plt\n'), ((93, 8, 93, 27), 'numpy.linspace', 'linspace', ({(93, 17, 93, 18): '0', (93, 20, 93, 21): '1', (93, 23, 93, 26): '501'}, {}), '(0, 1, 501)', False, 'from numpy import linspace\n'), ((99, 4, 99, 20), 'matplotlib.pyplot.plot', 'plt.plot', ({(99, 13, 99, 14): 'x', (99, 16, 99, 19): 'u_e'}, {}), '(x, u_e)', True, 'import matplotlib.pyplot as plt\n'), ((100, 4, 101, 33), 'matplotlib.pyplot.legend', 'plt.legend', (), '', True, 'import matplotlib.pyplot as plt\n'), ((102, 4, 102, 21), 'matplotlib.pyplot.xlabel', 'plt.xlabel', ({(102, 15, 102, 20): '"""$x$"""'}, {}), "('$x$')", True, 'import matplotlib.pyplot as plt\n'), ((102, 24, 102, 41), 'matplotlib.pyplot.ylabel', 'plt.ylabel', ({(102, 35, 102, 40): '"""$u$"""'}, {}), "('$u$')", True, 'import matplotlib.pyplot as plt\n'), ((104, 4, 104, 27), 'matplotlib.pyplot.savefig', 'plt.savefig', ({(104, 16, 104, 26): '"""tmpd.png"""'}, {}), "('tmpd.png')", True, 'import matplotlib.pyplot as plt\n'), ((104, 29, 104, 52), 'matplotlib.pyplot.savefig', 'plt.savefig', ({(104, 41, 104, 51): '"""tmpd.pdf"""'}, {}), "('tmpd.pdf')", True, 'import matplotlib.pyplot as plt\n'), ((118, 4, 118, 16), 'matplotlib.pyplot.figure', 'plt.figure', ({}, {}), '()', True, 'import matplotlib.pyplot as plt\n'), ((120, 4, 120, 14), 'matplotlib.pyplot.show', 'plt.show', ({}, {}), '()', True, 'import matplotlib.pyplot as plt\n'), ((7, 10, 7, 37), 'sympy.integrate', 'sym.integrate', ({(7, 24, 7, 25): '(1)', (7, 27, 7, 36): '(x, 0, x)'}, {}), '(1, (x, 0, x))', True, 'import sympy as sym\n'), ((8, 8, 8, 37), 'sympy.integrate', 'sym.integrate', ({(8, 22, 8, 25): 'u_x', (8, 27, 8, 36): '(x, 0, x)'}, {}), '(u_x, (x, 0, x))', True, 'import sympy as sym\n'), ((13, 21, 13, 34), 'sympy.expand', 'sym.expand', ({(13, 32, 13, 33): 'u'}, {}), '(u)', True, 'import sympy as sym\n'), ((48, 8, 48, 25), 'matplotlib.pyplot.plot', 'plt.plot', ({(48, 17, 48, 18): 'x', (48, 20, 48, 24): 'u[k]'}, {}), '(x, u[k])', True, 'import matplotlib.pyplot as plt\n'), ((98, 8, 98, 25), 'matplotlib.pyplot.plot', 'plt.plot', ({(98, 17, 98, 18): 'x', (98, 20, 98, 24): 'u[k]'}, {}), '(x, u[k])', True, 'import matplotlib.pyplot as plt\n'), ((37, 39, 37, 58), 'numpy.sin', 'sin', ({(37, 43, 37, 57): '((2 * i + 1) * pi * x / 2)'}, {}), '((2 * i + 1) * pi * x / 2)', False, 'from numpy import pi, sin, zeros\n'), ((83, 43, 83, 60), 'numpy.sin', 'sin', ({(83, 47, 83, 59): '((i + 1) * pi * x / 2)'}, {}), '((i + 1) * pi * x / 2)', False, 'from numpy import pi, sin, zeros\n'), ((28, 24, 28, 43), 'numpy.sin', 'sin', ({(28, 28, 28, 42): '(2 * i + 1) * pi * x / 2'}, {}), '((2 * i + 1) * pi * x / 2)', False, 'from numpy import pi, sin, zeros\n'), ((65, 24, 65, 41), 'numpy.sin', 'sin', ({(65, 28, 65, 40): '(i + 1) * pi * x / 2'}, {}), '((i + 1) * pi * x / 2)', False, 'from numpy import pi, sin, zeros\n'), ((73, 24, 73, 41), 'numpy.sin', 'sin', ({(73, 28, 73, 40): '(i + 1) * pi * x / 2'}, {}), '((i + 1) * pi * x / 2)', False, 'from numpy import pi, sin, zeros\n'), ((85, 43, 85, 60), 'numpy.sin', 'sin', ({(85, 47, 85, 59): '((i + 1) * pi * x / 2)'}, {}), '((i + 1) * pi * x / 2)', False, 'from numpy import pi, sin, zeros\n'), ((10, 19, 10, 32), 'sympy.diff', 'sym.diff', ({(10, 28, 10, 29): 'u', (10, 30, 10, 31): 'x'}, {}), '(u, x)', True, 'import sympy as sym\n')]