0
$\begingroup$

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]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.