library(ggplot2)
library(gganimate)
library(tidyr)

Problem 8.1

a)

We can approximate:

\[ \frac{\partial^2 u}{\partial t^2}=v^2\frac{\partial^2 u}{\partial x^2} \] By recalling (eqn. 8.20) that

\[ \frac{\partial^2u}{\partial x^2}=\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{(\Delta x)^2} \] and \[ \frac{\partial^2u}{\partial t^2}=\frac{u_j^{n+1}-2u_j^n+u_j^{n-1}}{(\Delta t)^2} \] Therefore… \[ \frac{u_j^{n+1}-2u_j^n+u_j^{n-1}}{(\Delta t)^2}=v^2\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{(\Delta x)^2} \] \[ u_j^{n+1}-2u_j^n+u_j^{n-1}=(\frac{v^2(\Delta t)^2}{(\Delta x)^2})(u_{j+1}^n-2u_j^n+u_{j-1}^n) \] \[ u_j^{n+1}=(\frac{v^2(\Delta t)^2}{(\Delta x)^2})(u_{j+1}^n-2u_j^n+u_{j-1}^n)+2u_j^n-u_j^{n-1} \]

b)

As one can see, this is a second order approximation in both time and space.

c)

Assume that \(u_n^j=A^ne^{ik(\Delta x)j}\) and substitute in to the above equation…

\[ A^{n+1}e^{ik(\Delta x)j}-2A^ne^{ik(\Delta x)j}+A^{n-1}e^{ik(\Delta x)j}=\frac{v^2(\Delta t)^2}{(\Delta x)^2}[A^ne^{ik(\Delta x)(j+1)}-A^ne^{ik(\Delta x)j}+A^ne^{ik(\Delta x)(j-1)}] \] Dividing through by \(A^n\) and \(e^{ik(\Delta x)j}\) and rearranging…

\[ A+\frac{1}{A}=\frac{v^2(\Delta t)^2}{(\Delta x)^2}(e^{ik\Delta x}-2+e^{-ik\Delta x})+2 \] \[ A+\frac{1}{A}=2+\frac{2v^2(\Delta t)^2}{(\Delta x)^2}(cos(k\Delta x)-1) \] As one can see, this is quadratic, as multiplying by \(A^2\) and rearranging gets us……

\[ A^2+A(2+\frac{2v^2(\Delta t)^2}{(\Delta x)^2}(cos(k\Delta x)-1))+1=0 \] so let’s use that to find A… Say that \(\alpha=2+\frac{2v^2(\Delta t)^2}{(\Delta x)^2}(cos(k\Delta x)-1)\)

\[ A=\frac{-\alpha\pm\sqrt{\alpha^2-4}}{2} \] \[ \to A=-1-\frac{v^2(\Delta t)^2}{(\Delta x)^2}\pm\sqrt{\frac{v^4(\Delta t)^4}{(\Delta x)^4}(cos(k\Delta x)-1)^2+2(cos(k\Delta x)-1)} \] The mode aplitudes are therefore \(A_1=\frac{-\alpha+\sqrt{\alpha^2-4}}{2}\) and \(A_2=\frac{-\alpha-\sqrt{\alpha^2-4}}{2}\)

d)

Following the hint, let us multiply together \(A_1\) and \(A_2\)…

\[ A_1A_2=\frac{-\alpha+\sqrt{\alpha^2-4}}{2}\frac{-\alpha-\sqrt{\alpha^2-4}}{2} \] \[ =\frac{\alpha^2}{4}+(\frac{1}{4}(\alpha^2-4))=1 \] Therfore \(A_1A_2=1\). What does this mean? It requires that they are complex conjugates with the term in the \(-\alpha/2\) as Real and \(\sqrt{\alpha^2-4}\) as imaginary, since that is the only way we could gaurantee that they are equal to 1. Thus, the square root term must be less than 0.

\[ \alpha^2-4<0 \] \[ \alpha^2<4 \] \[ -2<\alpha<2 \] Substituting in for \(\alpha\)…

\[ -2<2+\frac{2v^2(\Delta t)^2}{(\Delta x)^2}(cos(k\Delta x)-1)<2 \] \[ -4<\frac{2v^2(\Delta t)^2}{(\Delta x)^2}(cos(k\Delta x)-1)<0 \] \[ -2<\frac{v^2(\Delta t)^2}{(\Delta x)^2}(cos(k\Delta x)-1)<0 \] The parenthetical term must always be between \((-1,0)\). Therefore, the restriction on \(v,\Delta t,\) and \(\Delta x\) is that:

\[ \frac{v^2(\Delta t)^2}{(\Delta x)^2}<2 \] ## e)

As we can see, \(A_1\) and \(A_2\) have the same real parts, but different imaginary parts. Therefore, the real parts decay at the same rate, and thus the different modes decay at different rates.

#8.2

a)

D<-1
delta_x<-1
delta_t<-c(1,.5,.1)
X<-seq(1,500,delta_x)



diff_t<-list()

for(j in 1:length(delta_t)){
  u<-rep(0,length(X))
  u[length(u)/2]<-1
  results<-NULL
  results<-rbind(results,u)
  TT<-seq(0,100*pi,delta_t[j])
  
  for(t in 1:length(TT)){
    u<-results[t,]
    j_plus<-c(0,u)[1:500]
    j_plus[500]<-0
    j_minus<-c(u[2:500],0)
    j_minus[1]<-0
    
    temp<-u+(D*delta_t[j]/delta_x^2)*(j_plus-2*u+j_minus)
    results<-rbind(results,temp)
  }
  diff_t[[j]]<-results
}

knitr::kable(as.data.frame(diff_t[[1]])[1:10,240:260])
V240 V241 V242 V243 V244 V245 V246 V247 V248 V249 V250 V251 V252 V253 V254 V255 V256 V257 V258 V259 V260
u 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
temp 0 0 0 0 0 0 0 0 0 1 -1 1 0 0 0 0 0 0 0 0 0
temp.1 0 0 0 0 0 0 0 0 1 -2 3 -2 1 0 0 0 0 0 0 0 0
temp.2 0 0 0 0 0 0 0 1 -3 6 -7 6 -3 1 0 0 0 0 0 0 0
temp.3 0 0 0 0 0 0 1 -4 10 -16 19 -16 10 -4 1 0 0 0 0 0 0
temp.4 0 0 0 0 0 1 -5 15 -30 45 -51 45 -30 15 -5 1 0 0 0 0 0
temp.5 0 0 0 0 1 -6 21 -50 90 -126 141 -126 90 -50 21 -6 1 0 0 0 0
temp.6 0 0 0 1 -7 28 -77 161 -266 357 -393 357 -266 161 -77 28 -7 1 0 0 0
temp.7 0 0 1 -8 36 -112 266 -504 784 -1016 1107 -1016 784 -504 266 -112 36 -8 1 0 0
temp.8 0 1 -9 45 -156 414 -882 1554 -2304 2907 -3139 2907 -2304 1554 -882 414 -156 45 -9 1 0 
knitr::kable(as.data.frame(diff_t[[2]])[1:10,240:260])
V240 V241 V242 V243 V244 V245 V246 V247 V248 V249 V250 V251 V252 V253 V254 V255 V256 V257 V258 V259 V260
u 0 0.0000000 0.0000000 0.0000000 0.000000 0.0000000 0.000000 0.0000000 0.000000 0.0000000 1.0000000 0.0000000 0.000000 0.0000000 0.000000 0.0000000 0.000000 0.0000000 0.0000000 0.0000000 0
temp 0 0.0000000 0.0000000 0.0000000 0.000000 0.0000000 0.000000 0.0000000 0.000000 0.5000000 0.0000000 0.5000000 0.000000 0.0000000 0.000000 0.0000000 0.000000 0.0000000 0.0000000 0.0000000 0
temp.1 0 0.0000000 0.0000000 0.0000000 0.000000 0.0000000 0.000000 0.0000000 0.250000 0.0000000 0.5000000 0.0000000 0.250000 0.0000000 0.000000 0.0000000 0.000000 0.0000000 0.0000000 0.0000000 0
temp.2 0 0.0000000 0.0000000 0.0000000 0.000000 0.0000000 0.000000 0.1250000 0.000000 0.3750000 0.0000000 0.3750000 0.000000 0.1250000 0.000000 0.0000000 0.000000 0.0000000 0.0000000 0.0000000 0
temp.3 0 0.0000000 0.0000000 0.0000000 0.000000 0.0000000 0.062500 0.0000000 0.250000 0.0000000 0.3750000 0.0000000 0.250000 0.0000000 0.062500 0.0000000 0.000000 0.0000000 0.0000000 0.0000000 0
temp.4 0 0.0000000 0.0000000 0.0000000 0.000000 0.0312500 0.000000 0.1562500 0.000000 0.3125000 0.0000000 0.3125000 0.000000 0.1562500 0.000000 0.0312500 0.000000 0.0000000 0.0000000 0.0000000 0
temp.5 0 0.0000000 0.0000000 0.0000000 0.015625 0.0000000 0.093750 0.0000000 0.234375 0.0000000 0.3125000 0.0000000 0.234375 0.0000000 0.093750 0.0000000 0.015625 0.0000000 0.0000000 0.0000000 0
temp.6 0 0.0000000 0.0000000 0.0078125 0.000000 0.0546875 0.000000 0.1640625 0.000000 0.2734375 0.0000000 0.2734375 0.000000 0.1640625 0.000000 0.0546875 0.000000 0.0078125 0.0000000 0.0000000 0
temp.7 0 0.0000000 0.0039062 0.0000000 0.031250 0.0000000 0.109375 0.0000000 0.218750 0.0000000 0.2734375 0.0000000 0.218750 0.0000000 0.109375 0.0000000 0.031250 0.0000000 0.0039062 0.0000000 0
temp.8 0 0.0019531 0.0000000 0.0175781 0.000000 0.0703125 0.000000 0.1640625 0.000000 0.2460938 0.0000000 0.2460938 0.000000 0.1640625 0.000000 0.0703125 0.000000 0.0175781 0.0000000 0.0019531 0 
knitr::kable(as.data.frame(diff_t[[3]])[1:10,240:260])
V240 V241 V242 V243 V244 V245 V246 V247 V248 V249 V250 V251 V252 V253 V254 V255 V256 V257 V258 V259 V260
u 0 0 0e+00 0.0e+00 0.00e+00 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00e+00 0.0e+00 0e+00 0 0
temp 0 0 0e+00 0.0e+00 0.00e+00 0.0000000 0.0000000 0.0000000 0.0000000 0.1000000 0.8000000 0.1000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00e+00 0.0e+00 0e+00 0 0
temp.1 0 0 0e+00 0.0e+00 0.00e+00 0.0000000 0.0000000 0.0000000 0.0100000 0.1600000 0.6600000 0.1600000 0.0100000 0.0000000 0.0000000 0.0000000 0.00e+00 0.0e+00 0e+00 0 0
temp.2 0 0 0e+00 0.0e+00 0.00e+00 0.0000000 0.0000000 0.0010000 0.0240000 0.1950000 0.5600000 0.1950000 0.0240000 0.0010000 0.0000000 0.0000000 0.00e+00 0.0e+00 0e+00 0 0
temp.3 0 0 0e+00 0.0e+00 0.00e+00 0.0000000 0.0001000 0.0032000 0.0388000 0.2144000 0.4870000 0.2144000 0.0388000 0.0032000 0.0001000 0.0000000 0.00e+00 0.0e+00 0e+00 0 0
temp.4 0 0 0e+00 0.0e+00 0.00e+00 0.0000100 0.0004000 0.0064500 0.0528000 0.2241000 0.4324800 0.2241000 0.0528000 0.0064500 0.0004000 0.0000100 0.00e+00 0.0e+00 0e+00 0 0
temp.5 0 0 0e+00 0.0e+00 1.00e-06 0.0000480 0.0009660 0.0104800 0.0652950 0.2278080 0.3908040 0.2278080 0.0652950 0.0104800 0.0009660 0.0000480 1.00e-06 0.0e+00 0e+00 0 0
temp.6 0 0 0e+00 1.0e-07 5.60e-06 0.0001351 0.0018256 0.0150101 0.0760648 0.2278563 0.3582048 0.2278563 0.0760648 0.0150101 0.0018256 0.0001351 5.60e-06 1.0e-07 0e+00 0 0
temp.7 0 0 0e+00 6.0e-07 1.80e-05 0.0002912 0.0029750 0.0197971 0.0851385 0.2257120 0.3321351 0.2257120 0.0851385 0.0197971 0.0029750 0.0002912 1.80e-05 6.0e-07 0e+00 0 0
temp.8 0 0 1e-07 2.3e-06 4.36e-05 0.0005323 0.0043888 0.0246490 0.0926617 0.2222970 0.3108505 0.2222970 0.0926617 0.0246490 0.0043888 0.0005323 4.36e-05 2.3e-06 1e-07 0 0

I will animate and finish problem set next week…