# How to implement the following Finite Element method for Burgers' equation?

I am trying to replicate this result.

It involves using the Galerkin finite element approach onto the viscous Burgers' equation. However, my implementation (in R) seems to be giving me wrong results. For example, if you refer to Table 1, $$u(x=0.25, t=0.25)=0.550$$ but my calculations give me $$0.001318505$$. It is completely off the mark.

I know some of my mistakes are due to incorrectly implementing $$u[m,n]$$ and not utilizing the boundary conditions correctly.

In my implementation, my main confusion is how to obtain the approximation $$u[m,n]$$, where $$m$$ represents the spatial points and $$n$$ represents the time $$t$$. If I am understood the article correctly, the initial condition is given as

$$U(x,0)=\sin(\pi x), x\in[0,1]$$

but I don't understand how to relate $$U$$ back to $$u$$. So what I did was just treat $$U=u$$.

Secondly, the boundary conditions are given as

$$U(x=0,t)=U(x=1,t)=0, t>0$$

So what I did was

1. When $$u[m==1,n]$$, I removed all $$u[m-1,n]$$, i.e. setting them to $$0$$

2. When $$u[m==N+1,n]$$, I removed all $$u[m+1,n]$$, i.e. setting them to $$0$$

So how should I correct my mistakes?

My code in R

v=1
# h=c(0.25, 0.125, 0.0625)
h=0.25
b=v/h[1]^2

input_dt=0.01
T=0.25/input_dt

x_values=seq(from=0, to=1, by=h[1])
N=(1/h[1])

A[i,j]=(1/6)*(matrix(c(2,1,1,2), ncol=2,nrow=2))
C[i,j]=1/2*(u_1/h)*matrix(c(-1,-1,1,1), ncol=2,nrow=2)
D[i,j]=matrix(c(1,-1,-1,1), ncol=2,nrow=2)

u=matrix(0, ncol=T+1, nrow=N+1)
u[,1]=sin(pi*x_values)

for(n in 1:(T)){
A=numeric()
B=numeric()
C=numeric()
D=numeric()
for(m in 1:(N+1)){
if(m==1){
# Remove u[m-1,] (as 0) when m==1
# LHS
A[m]=1/6 - b*input_dt/2 - (input_dt/4)*(0) # Remove u[m-1,] (as 0) when m==1
B[m]=2/3 + b*input_dt + (input_dt/4)*((0-u[m,n])/h)
C[m]=1/6 - b*input_dt/2 + (input_dt/4)*(u[m,n]/h)
#RHS
D[m]=(1/6 - b*input_dt/2 - (input_dt/4)*((0)/h))*0 + (2/3-b*input_dt-(input_dt/4)*(0/h-u[m,n]/h))*u[m,n]+(1/6+b*input_dt/2-(input_dt/4)*u[m,n])*u[m+1,n]
}else if(m==(N+1)){
# Remove u[m+1,] (as 0) when m==N+1
# LHS
A[m]=1/6 - b*input_dt/2 - (input_dt/4)*(u[m-1,n]/h)
B[m]=2/3 + b*input_dt + (input_dt/4)*((u[m-1,n]-u[m,n])/h)
C[m]=1/6 - b*input_dt/2 + (input_dt/4)*(u[m,n]/h)
#RHS
D[m]=(1/6 - b*input_dt/2 - (input_dt/4)*((u[m-1,n])/h))*u[m-1,n] + (2/3-b*input_dt-(input_dt/4)*(u[m-1,n]/h-u[m,n]/h))*u[m,n]+(1/6+b*input_dt/2-(input_dt/4)*u[m,n])*0
}else{
# LHS
A[m]=1/6 - b*input_dt/2 - (input_dt/4)*(u[m-1,n]/h)
B[m]=2/3 + b*input_dt + (input_dt/4)*((u[m-1,n]-u[m,n])/h)
C[m]=1/6-(b*input_dt/2) + (input_dt/4)*(u[m,n]/h)
#RHS
D[m]=(1/6 - b*input_dt/2 - (input_dt/4)*((u[m-1,n])/h))*u[m-1,n] + (2/3-b*input_dt-(input_dt/4)*(u[m-1,n]/h-u[m,n]/h))*u[m,n]+(1/6+b*input_dt/2-(input_dt/4)*u[m,n])*u[m+1,n]
}
}

# Set up tridiagonal-matrix
mat = diag(B)
mat[row(mat) - col(mat) == 1]  = A[-1]
mat[row(mat) - col(mat) == -1] = C[-length(C)]

# Solve for solutions
sol=(base::solve(mat))%*%t(t(D))

# Update solution matrix
u[ ,(n+1)]=sol

}



To get $0.001318505$

u[2,26]