Plotting a solution of a heat equation IBVP

In [52]:
%pylab inline
x = linspace(0,2*pi,11)
y = sin(x)
plot(x,y);

Populating the interactive namespace from numpy and matplotlib

We draw smooth-looking curves by drawing polygonal curves with very many points:

In [51]:
x = linspace(0,2*pi,211)
y = sin(x)
plot(x,y);
In [5]:
x = linspace(0,2*pi,5)
x
Out[5]:
array([ 0.        ,  1.57079633,  3.14159265,  4.71238898,  6.28318531])
In [6]:
100*x
Out[6]:
array([   0.        ,  157.07963268,  314.15926536,  471.23889804,
        628.31853072])
In [7]:
y = sin(x)
y
Out[7]:
array([  0.00000000e+00,   1.00000000e+00,   1.22464680e-16,
        -1.00000000e+00,  -2.44929360e-16])
In [9]:
L = 2.
t = .01
n = 12
x = linspace(0,L,200)
y = sin(n*pi*x/L)*exp(-(n*pi/L)**2*t )
plot(x,y)
Out[9]:
[<matplotlib.lines.Line2D at 0x7f6c7afe9c18>]
In [10]:
for t in [1,2,5]:
    print(t)

1
2
5
In [12]:
for t in linspace(0,.1,6):
    print(t)

0.0
0.02
0.04
0.06
0.08
0.1

First plot just one of our elementary solutions:

In [15]:
n = 12
for t in linspace(0,.01,6):
    y = sin(n*pi*x/L)*exp(-(n*pi/L)**2*t )
    plot(x,y)
    print(t)

0.0
0.002
0.004
0.006
0.008
0.01

Then another one, with a smaller value of n:

In [16]:
n = 4
for t in linspace(0,.01,6):
    y = sin(n*pi*x/L)*exp(-(n*pi/L)**2*t )
    plot(x,y)
    print(t)

0.0
0.002
0.004
0.006
0.008
0.01

A linear combination of two of the elementary solutions:

In [19]:
for t in linspace(0,.01,6):
    y = 1.3*sin(4*pi*x/L)*exp(-(4*pi/L)**2*t ) + 0.5*sin(12*pi*x/L)*exp(-(12*pi/L)**2*t )
    plot(x,y)
    print(t)

0.0
0.002
0.004
0.006
0.008
0.01

Now we will build the solution of the IBVP we figured out in class last time.

In [22]:
M = 4
x = linspace(0,L,5)
n = arange(1,M+1)
n
Out[22]:
array([1, 2, 3, 4])

We will build a table of values of $\sin\frac{n\pi x}{L}$, in stages ...

In [23]:
outer(x,n)
Out[23]:
array([[ 0. ,  0. ,  0. ,  0. ],
       [ 0.5,  1. ,  1.5,  2. ],
       [ 1. ,  2. ,  3. ,  4. ],
       [ 1.5,  3. ,  4.5,  6. ],
       [ 2. ,  4. ,  6. ,  8. ]])

If matrix o = outer(x,n), then $$o_{ij} = x_i n_j$$

In [24]:
10*outer(x,n)
Out[24]:
array([[  0.,   0.,   0.,   0.],
       [  5.,  10.,  15.,  20.],
       [ 10.,  20.,  30.,  40.],
       [ 15.,  30.,  45.,  60.],
       [ 20.,  40.,  60.,  80.]])
In [25]:
outer(x,n)*pi/L
Out[25]:
array([[  0.        ,   0.        ,   0.        ,   0.        ],
       [  0.78539816,   1.57079633,   2.35619449,   3.14159265],
       [  1.57079633,   3.14159265,   4.71238898,   6.28318531],
       [  2.35619449,   4.71238898,   7.06858347,   9.42477796],
       [  3.14159265,   6.28318531,   9.42477796,  12.56637061]])
In [27]:
sin(outer(x,n)*pi/L).round(5)
Out[27]:
array([[ 0.     ,  0.     ,  0.     ,  0.     ],
       [ 0.70711,  1.     ,  0.70711,  0.     ],
       [ 1.     ,  0.     , -1.     , -0.     ],
       [ 0.70711, -1.     ,  0.70711,  0.     ],
       [ 0.     , -0.     ,  0.     , -0.     ]])
In [28]:
B = 200*(1-cos(n*pi))/(n*pi)
B
Out[28]:
array([ 127.32395447,    0.        ,   42.44131816,    0.        ])
In [30]:
(B*sin(outer(x,n)*pi/L)).round(5)
Out[30]:
array([[   0.     ,    0.     ,    0.     ,    0.     ],
       [  90.03163,    0.     ,   30.01054,    0.     ],
       [ 127.32395,    0.     ,  -42.44132,   -0.     ],
       [  90.03163,   -0.     ,   30.01054,    0.     ],
       [   0.     ,   -0.     ,    0.     ,   -0.     ]])
In [33]:
for t in linspace(0,.01,6):
    Q = (B*sin(outer(x,n)*pi/L))*exp(-(n*pi/L)**2*t)
    u = Q.sum(axis=1)
    print(u.round(5))

[   0.       120.04218   84.88264  120.04218    0.     ]
[   0.       118.29528   86.09958  118.29528    0.     ]
[   0.       116.6072    87.23951  116.6072     0.     ]
[   0.       114.97547   88.30589  114.97547    0.     ]
[   0.       113.39772   89.30205  113.39772    0.     ]
[   0.       111.87169   90.23115  111.87169    0.     ]
In [34]:
for t in linspace(0,.01,6):
    Q = (B*sin(outer(x,n)*pi/L))*exp(-(n*pi/L)**2*t)
    u = Q.sum(axis=1)
    plot(x,u)

Now increase the number of x values to make the curves smooth, and increase the number of included terms, M, to make the approximations more accurate:

In [48]:
M = 40
tfinal = 1.
x = linspace(0,L,200)
n = arange(1,M+1)
B = 200*(1-cos(n*pi))/(n*pi)
for t in linspace(0,tfinal,100)[1:]:
    Q = (B*sin(outer(x,n)*pi/L))*exp(-(n*pi/L)**2*t)
    u = Q.sum(axis=1)
    plot(x,u,'b',alpha=0.5)
xlabel('x')
ylabel('u(x,t)');
In [ ]: