(9.1) Consider the damped wave equation
$$\frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2} + \gamma \frac{\partial}{\partial t} \frac{\partial^2 u}{\partial x^2}
$$
Take the solution domain to be the interval [0,1].
(a) Use the Galerkin method to find an approximating system of differenctial equations.
(b) Evaluate the matrix coefficent for linear hat functions, using elements with a fixed size of h.
(c) Now find the matrix coefficient for Hermit polynomial interpolation basis functions, once again using elements with a fixed size of h. A symbolic math enviroment is useful for this problem.
#%matplotlib inline
#import matplotlib.pyplot as plt
from sympy import *
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
#init_printing()
# or with older versions of sympy/ipython, load the IPython extension
#%load_ext sympy.interactive.ipythonprinting
# or
#%load_ext sympyprinting
Define Matrix
u_0,du_0,u_h,du_h,a_0,a_1,a_2,a_3,h,x,f,df,temp = symbols("u_0, du_0,u_h,du_h,a_0,a_1,a_2,a_3,h,x,f,df,temp ")
U = Matrix([[u_0],[du_0],[u_h],[du_h]])
M = Matrix([[1,0,0,0], [0,1,0,0],[1,h,h**2,h**3],[0,1,2*h,3*h**2]])
A = Matrix([[a_0],[a_1],[a_2],[a_3]])
X = Matrix([[1, x, x**2, x**3]])
X
f =[]
df =[]
Minv = M.inv()
##print(Minv)
for i in range(0,4):
temp = 0
for j in range(0,4):
temp += Minv[:, i][j]*X[j]
f.append(temp)
df.append(diff(temp,x))
print(f)
print(df)
For coeff A and B:
1) when |i-j| = 0, $\phi_i = u_3,\phi_j = u_3 $
2) when |i-j| = 0, $\phi_i = u_3,\phi_j = u_4 $
3) when |i-j| = 0, $\phi_i = u_4,\phi_j = u_4 $
4) when |i-j| = 1, $\phi_i = u_3,\phi_j = u_3 $
5) when |i-j| = 1, $\phi_i = u_3,\phi_j = u_4 $
6) when |i-j| = 1, $\phi_i = u_4,\phi_j = u_4 $
7) when |i-j| > 1
## integrate(f[2]*f[2], (x, 0, 2*h)) // if we do so, we need to define M value again , since boundary condition is changed. So, we use rather shifted graph for get integration.
A_0_3_3 = integrate(f[2]*f[2], (x, 0, h)) + integrate(f[0]*f[0], (x, 0, h))
A_0_3_4 = integrate(f[2]*f[3], (x, 0, h)) + integrate(f[0]*f[1], (x, 0, h))
A_0_4_4 = integrate(f[3]*f[3], (x, 0, h)) + integrate(f[1]*f[1], (x, 0, h))
A_1_3_3 = integrate(f[0]*f[2], (x, 0, h))
A_1_3_4 = integrate(f[1]*f[2], (x, 0, h))
A_1_4_4 = integrate(f[1]*f[3], (x, 0, h))
B_0_3_3 = integrate(df[2]*df[2], (x, 0, h)) + integrate(df[0]*df[0], (x, 0, h))
B_0_3_4 = integrate(df[2]*df[3], (x, 0, h)) + integrate(df[0]*df[1], (x, 0, h))
B_0_4_4 = integrate(df[3]*df[3], (x, 0, h)) + integrate(df[1]*df[1], (x, 0, h))
B_1_3_3 = integrate(df[0]*df[2], (x, 0, h))
B_1_3_4 = integrate(df[1]*df[2], (x, 0, h))
B_1_4_4 = integrate(df[1]*df[3], (x, 0, h))
print("Matrix_|i-j|_fuction1_fuction2 : "+ "coefficient value")
print("A_0_3_3 : "+ str(A_0_3_3))
print("A_0_3_4 : "+ str(A_0_3_4))
print("A_0_4_4 : "+ str(A_0_4_4))
print("A_1_3_3 : "+ str(A_1_3_3))
print("A_1_3_4 : "+ str(A_1_3_4))
print("A_1_4_4 : "+ str(A_1_4_4))
print("B_0_3_3 : "+ str(B_0_3_3))
print("B_0_3_4 : "+ str(B_0_3_4))
print("B_0_4_4 : "+ str(B_0_4_4))
print("B_1_3_3 : "+ str(B_1_3_3))
print("B_1_3_4 : "+ str(B_1_3_4))
print("B_1_4_4 : "+ str(B_1_4_4))
Now Let's implement this to wave equation
## let's use both 3,3
## if |i-j|=0, value is Aij = 26*h/35, Bij = 12/(5*h)
## if |i-j|=1, value is Aij = 9*h/70, Bij = -6/(5*h)
h = 0.1
boundary = 1
count = int(boundary/h)
A = []
B = []
for j in range (0,count):
tempA = []
tempB = []
for i in range(0,count):
if (abs(i-j)==0 ):
##tempA.append(26*h/35)
##tempB.append(12/(5*h))
tempA.append(2*h/3)
tempB.append(2/h)
elif (abs(i-j)==1):
##tempA.append(9*h/70)
##tempB.append(-6/(5*h))
tempA.append(1*h/6)
tempB.append(-1/h)
else:
tempA.append(0)
tempB.append(0)
A.append(tempA)
B.append(tempB)
M_A = Matrix(A)
M_B = Matrix(B)
##M_A_inv = M_A.inv()
##Result = M_B*M_A_inv
##Result
From 9.1 (a),
$$ A_{ij}\cdot \frac{d^2 \vec{a}}{d t^2} + \gamma\cdot B_{ij} \cdot \frac{d\vec{a}}{dt} + v^2\cdot B_{ij} \cdot \vec{a} = \vec{0}$$when at certain t, from equation(9.9)
$$u(x_{i}) = a_i $$
so vector a is value of u in certain time. $$\vec{a} = \vec{u}$$
From partial difference chapter(8) we got,
$$ \frac{\partial u}{\partial t} = \frac{u_{j}^{n+1} - u_j^n}{{\Delta t}} $$
$$ \frac{\partial^2 u}{\partial t^2} = \frac{u_{j}^{n+1} - 2 u_j^n + u_{j}^{n-1}}{({\Delta t})^2}$$
So, 9.1(a) is again,
$$A \cdot (\vec{u}^{n+1}-2\vec{u}^{n}+\vec{u}^{n-1}) + \gamma B \cdot dt \cdot (\vec{u}^{n+1} - \vec{u}^{n}) + v^2 B \cdot dt^2 \cdot \vec{u}^n = \vec{0}$$So, this indicates, if we know $u^0$ and $u^1$, we can get $u^n$.
Let's
$$\vec{u}^{n+1} =M_1[M_2\cdot \vec{u}^n - A\cdot \vec{u}^{n-1}]$$
## let's make a wave
dt = 0.01
gamma = 0
v = 1
##h = 1
M_1 = (M_A + M_B*gamma*dt).inv()
M_2 = 2*M_A + gamma*dt*M_B -(v**2)*(dt**2)*M_B
M_1
M_2
u_0=[]
u_1=[]
for i in range (0,count):
temp_u_0 = [0]
temp_u_1 = [0]
u_0.append(temp_u_0)
u_1.append(temp_u_1)
M_u_0 = Matrix(u_0)
M_u_1 = Matrix(u_1)
## update one value
M_u_1[5,0] = 5
M_u_1
M_u_2 = M_1*(M_2*M_u_1 - M_A*M_u_0)
M_u_0 = M_u_1
M_u_1 = M_u_2
##current
current = M_u_1.T.tolist()[0]
x=[]
for i in range(0, count):
x.append(i)
plt.plot(x, current,'b',c='red',label="theta")
plt.xlabel("x")
plt.ylabel("y")
plt.show()
