1
$\begingroup$

I am attempting to replicate results from this article.

I'm not sure why but my results are completely different and wrong. For example, the exact solution with parameters ($x=0.1$, $T=0.01$, $Re=0.1$ and $\Delta t=0.0001$) is $0.11461$ but my approximated solution gives $-175.4007$.

Is anyone able to spot my mistake? I've included my code (in R) below.

The Thomas algorithm was implemented.

I've asked a very similar question in the Math SE website recently but I think with the focus on the coding as well, this perhaps might be a better place to post this question.


In case you can't access the article, these are the relevant information:

$0<l<N$ represents spatial point, $0<j<M$ represents time, $k=T/M$, $h=1/N$ and $x_i=ih$.

The initial condition for $0<x<1$ is given as $$\phi(x,0)=\exp\left({-\frac{Re}{2\pi}\left[1-\cos(\pi x)\right]}\right)$$

and boundary conditions for $0\leq t \leq T$ are $$\phi_x(0,t)=\phi_x(1,t)=0.$$

Then the Crank-Nicolson finite difference discretization equations are given as:

For $l=0$ $$-r\phi_{l+1,j+1}+(1+r)\phi_{l,j+1}+0=r\phi_{l+1,j}+(1-r)\phi_{l,j}$$

For $l=1:N-1$ $$-\frac{r}{2}\phi_{l+1,j+1}+(1+r)\phi_{l,j+1}-\frac{r}{2}\phi_{l-1,j+1}=\frac{r}{2}\phi_{l+1,j}+(1-r)\phi_{l,j}+\frac{r}{2}\phi_{l-1,j}$$

For $l=N$ $$0+(1+r)\phi_{l,j+1}-r\phi_{l-1,j+1}=r\phi_{l-1,j}+(1-r)\phi_{l,j}$$

where $r=k/Reh^2$.

The solution via Hopf-Cole transformation is then calculated as:

$$u_{l,j}=-\left(\frac{1}{Re}\right){\frac{\phi_{l+1,j}-\phi_{l-1,j}}{h\phi_{l,j}}}$$

# Parameters
Re=0.1
T=0.01
input_dt=0.0001
N=80

M=T/input_dt
h=1/N
k=T/M
r=k/(Re*(h^2))


# Spatial discretization
x_discrete=numeric()
x_discrete[1]=0
for( i in 1:(N+1)){
    x_discrete[i] = (i-1)*h
}
x_discrete[N+1]=1

# Calculate initial conditions
u_init=exp(-(Re/(2*pi))*(1-cos(pi*x_discrete)))
# Adjust u_init with BC
u_init[1]=0
u_init[N+1]=0

# matrix to store solutions
u=matrix(0, (N+1), (M+1))
# Add in initial conditions
u[,1]=u_init

# Crank-Nicolson
# For each time-step
for(j in 1:(M)){
    # A,B,C, D set up from:
    # https://www.sfu.ca/~rjones/bus864/notes/notes2.pdf
    A=numeric()
    B=numeric()
    C=numeric()
    D=numeric()

    # For each spatial node (x values)
    for(l in 1:(N+1)){
        if(l==1){
            # LHS
            A[l]=-r
            B[l]=(1+r)
            C[l]=0
            #RHS
            D[l]=r*u[l+1,j]+(1-r)*u[l,j]
        }else if(l==(N+1)){
            # LHS
            A[l]=0
            B[l]=(1+r)
            C[l]=-r
            #RHS
            D[l]=r*u[l-1,j]+(1-r)*u[l,j]
        }else{
            # LHS
            A[l]=-r/2
            B[l]=(1+r)
            C[l]=-r/2
            #RHS
            D[l]=(r/2)*u[l+1,j]+(1-r)*u[l,j]+(r/2)*u[l-1,j]
        }
    }

    # Set up tridiagonal-matrix
    # Method taken from:
    # https://stackoverflow.com/questions/49536400/create-a-tridiagonal-matrix
    m = diag(B)
    m[row(m) - col(m) == 1]  = A[-1]
    m[row(m) - col(m) == -1] = C[-length(C)]

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

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

    u[1, (j+1)]=0       # Adjust according to BC
    u[(N+1), (j+1)]=0       # Adjust according to BC
}

# Exact solution at x=0.1 should be 0.11461 (from Table 1)
# BUT mine gives -175 (from the line of code below)
-(1/Re)*(u[9+1,101] - u[9-1,101])/(h*u[9,101])
$\endgroup$
  • 1
    $\begingroup$ It would be helpful if you could add how the results are wrong, are they merely unexpectedly imprecise, or totally implausible. $\endgroup$ – Dr. Lutz Lehmann Aug 18 at 14:31
  • $\begingroup$ What exactly are your boundary conditions? In the end you correct the values to zero, but in the matrix construction you use that the derivatives are zero to set $u^{n}_{-1}=u^n_1$. This is contradicting and might be responsible for the error you observe. $\endgroup$ – Dr. Lutz Lehmann Aug 18 at 15:15
  • $\begingroup$ Thanks for taking the time to look at this. If possible, could you perhaps let me know which lines of code you are referring to? Did you mean the two lines at the bottom u[1,(j+1)]=0 and u[(N+1),(j+1)]=0 ? Is this correct? Since the boundary conditions state $\phi_x(x=0,t) = \phi_x(x=1,t)=0$, so all values at the start and end of each column ("spatial step") will be set to 0. Or am I not understanding the boundary conditions correctly? $\endgroup$ – mathnoob Aug 18 at 19:55
  • $\begingroup$ Yes, these lines. I read the BC as the $x$ derivatives being zero, so that you can continue a solution by its mirror image for one or two steps. $\endgroup$ – Dr. Lutz Lehmann Aug 18 at 22:08

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.