-1
$\begingroup$

I need to implement a numerical scheme for the solution of the one-dimensional advection equation $$\\\frac{\partial u}{\partial t} + C(x, t) \frac{\partial u}{\partial x} = 0 \\\\$$ $$ \\ C(x,t) = \frac{\pi\cos(2 \pi t) + 3.5}{2\pi - \pi\sin(\pi x)} \\\\$$ $$\\ 0 \leq x \leq 1, \quad 0 \leq t \leq 1\\\\$$

with Crank–Nicolson finite difference method $$\\ \frac{u^{j + 1}_{k} - u^{j}_{k}}{\tau} + \frac{C}{4h}(u^{j + 1}_{k + 1} - u^{j + 1}_{k - 1} + u^{j}_{k + 1} - u^{j}_{k - 1}) = 0\\\\$$ and boundary conditions $$\\ \phi(x) = u(x,0), \quad \psi_{0}(t) = u(0,t), \quad \psi_{1}(t) = u(1, t) \\\\$$

it's known that the exact solution has the following form $$\\ U(x,t) = \cos(\pi x) - \frac{1}{2}\sin(2 \pi t) + 2\pi x - 3.5t\\\\$$

The numerical solution on each layer in time is found as the solution of a system of linear equations with a tridiagonal matrix $$\\ Au = b \\\\$$ $$\\ A = \begin{pmatrix} 1 & \frac{C\tau}{4h} & 0 & ...\\ -\frac{C\tau}{4h} & 1 & \frac{C\tau}{4h} & ... \\ 0 & -\frac{C\tau}{4h} & 1 & ... \\ ... & ... & ... & ... \end{pmatrix}\\\\$$

I've write a simple python script, but the result doesn't match the theory:

graph

I'm sure I'm wrong somewhere, but I don't see where. Thanks!

$\endgroup$
1
  • 1
    $\begingroup$ Speed $C$ is positive at $x=1$. You cannot specify dirichlet conditions at $x=1$. $\endgroup$ – cfdlab Jan 5 at 13:47
0
$\begingroup$

Some remarks about the code, perhaps cleaning these up already solves the question, especially the initialization error:

Floating point errors are a thing, so if you want to get the expected results, give some wiggle room (0.01 is arbitrary, I would not expect defects above 1e-12 in this place)

  N = int((x1 - x0) / h + 0.01)
  M = int((t1 - t0) / tau + 0.01)

There is no guarantee that $N·h=x_1-x_0$ etc., so force this grid condition

h = (x1 - x0) / N
tau = (t1 - t0) / M

There is a honest error in the side boundary initialization, the space component is always first

    u[0, i] = psi_0(i * tau)
    u[N, i] = psi_1(i * tau)

The implicit trapezoidal method that is the method in time direction of C-N reads as $$ \frac{y^{j+1}-y^j}{\tau}=\frac{F(t^j,y^j)+F(t^{j+1},y^{j+1})}2. $$ Accordingly, the coefficients a_h for the banded matrix should be computed at time (j+1)*tau to get the correct order $O(\tau^2+h^2)$ of the method.

$\endgroup$
1
  • $\begingroup$ Thank you for the remarks, I am little ashamed of such obvious errors. Unfortunately, fixing these doesn't solve the question, perhaps there are more issues. $\endgroup$ – Winabryr Jan 5 at 19:40

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.