1
$\begingroup$

I would like to solve numerically the diffusion equation, where the sink term depends linearly on the field, and there is field-independent sink:

$\frac{\partial^2 u(x)}{\partial x^2} =f(x)u(x) - \kappa(x)$

At the moment I am solving the above using a finite difference scheme, where I modify the standard laplacian matrix, to encompass the $f(x)u(x)$ term in the diagonal components. I then solve a linear system of the form:

$Ax = b$

Where b is a vector composed of the discretised $\kappa(x)$ term, and $A$ is the standard finite difference laplacian (forgetting about boundary conditions), with modified diagonal components.

The above works fine, but the issue is that I am solving it for the quasi-steady-state at each time step (basically $f(x)$ and $\kappa(x)$ change over time). Since I am required to modify $A$ according to the linear source term at each time step, this means I am required to solve the system afresh each time.

What I would prefer is a method where I modify $b$ each time step, and keep $A$ constant. This would allow me to take advantage of $A$'s constancy, and do the bulk of solving once, and then just apply it at each time step.

I have tried a scheme whereby I use the current time step's value of $u_i$ on the RHS opposed to the future value, leading to a scheme of the form:

$\frac{u^{t+1}_{i+1} + u^{t+1}_{i-1} - 2u^{t+1}_i}{2h} = f_i u^t_i - \kappa_i$

However, this does not seem to work, and leads to instability/non-convergent behaviour.

Does anyone know of any schemes which would allow me to keep $A$ constant, and only modify $b$?

Alternatively, more generally can anyone tell me if there are other methods which might be better than finite difference in this case?

Best,

Ben

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ There is no time derivative in your original PDE, so why did you introduce a time-stepping scheme? $\endgroup$
    – Bill Barth
    Dec 4, 2014 at 13:18
  • $\begingroup$ @BillBarth - ok, the functions $f(x)$ and $\kappa(x)$ change over time, meaning that I need to resolve the above equation at specific points through time. Does that make sense? $\endgroup$
    – ben18785
    Dec 4, 2014 at 14:02
  • 1
    $\begingroup$ No, I don't see a $t$ in your PDE at all. Are you saying that your equation is actually $\frac{\partial^2 u(x,t)}{\partial x^2} = f(x,t)u(x,t) + \kappa(x,t)$? $\endgroup$
    – Bill Barth
    Dec 4, 2014 at 14:16
  • $\begingroup$ @BillBarth. Yes. $\endgroup$
    – ben18785
    Dec 4, 2014 at 15:15

3 Answers 3

Reset to default
2
$\begingroup$

In the comments you mentioned that $f$ has a very special form. Let $u:[a,b]\rightarrow \mathbb{R}$. Then $f:[a,b]\rightarrow\{0,1\}$ is defined for some set $S \subset [a,b]$ as follows: $$f(x) = \begin{cases} 1 & \text{if $x \in S$,} \\ 0 & \text{otherwise.} \end{cases}$$ As Bill Barth mentioned, you can now split the matrix $A$ into two parts as $$A = K+C.$$ I like to think in terms of finite element methods since I am more acquainted with those. Now your weak formulation is: find $u \in H^1_0(a,b)$ such that $$(u_x,v_x)+(fu,v)=(\kappa,v)$$ hold for very test function $v \in H^1_0(a,b)$.

Discretization of the first term leads to $K$ and the second term to $C$. My intuition now says that since your $f$ is so "simple", you can just build the mass matrix $M$ corresponding to the term $(u,v)$ and then construct a mapping $F_f$ satisfying $$C = F_f(M).$$ If $F_f$ is simple enough to evaluate (i.e. just setting rows and colums to zero), then you only once construct $M$ and every $C$ is just a slight perturbation.

Does it make any sense?

Improve this answer
$\endgroup$
1
  • $\begingroup$ thanks! I think I understand what you mean at a basic level. However, ultimately, I will still have to solve for a steady state by solving $Ax=b$, with A still being modified at each time step. Is there reason to believe that these modifications to $A$, and the resulting solving method are much faster than that which I am currently achieving with finite differences? $\endgroup$
    – ben18785
    Dec 8, 2014 at 16:37
1
$\begingroup$

You can't keep $A$ constant, but you could split it into two parts $A(t)=K+C(t)$ for the diffusion and non-constant "reaction" term and then use the $LU$-factorization of $K$ as a preconditioner for an iterative method for the overall system at future times. If $f(x,t)$ is reasonably nice, that should work well. Also, since there's no time-derivative in the PDE, each time level is independent, and you can solve many of them in parallel.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks, I think that may work. The only issue is that $f(x,t)$ is really not nice - it is really a discrete function where cell positions determine whether ($f(x,t)=1$,$\kappa=0$) or ($f(x,t)=0$,$\kappa=1$) or finally ($f(x,t)=0$,$\kappa=0$). Would this method still be applicable? $\endgroup$
    – ben18785
    Dec 4, 2014 at 15:44
  • $\begingroup$ That's not the best, but it's probably fine. It's definitely worth a shot. This all may be moot in 1D. The cost of a direct solve on a tri-dagonal matrix should be very competitive with even a well-preconditioned Krylov solver. I think you're better off parallelizing in time. $\endgroup$
    – Bill Barth
    Dec 4, 2014 at 15:50
  • $\begingroup$ Thank you. I am not sure if I can parallelize in time though, since $f$ and $\kappa$ are essentially Markov processes depending on the past time period? I will have a look at well-conditioned Krylov solvers though. $\endgroup$
    – ben18785
    Dec 4, 2014 at 16:19
  • $\begingroup$ Maybe you should put that part about $f$ and $\kappa$ in the question. You keep leaving out important details. $\endgroup$
    – Bill Barth
    Dec 4, 2014 at 16:22
1
$\begingroup$

You are basically trying to solve a one-dimensional modified Helmholtz equation with variable coefficients.

If the locations where the changes in the coefficients occur remain constant in time, I would suggest using a standard Finite Element (FE) approach where you would place element interfaces at these locations. Piecewise constant coefficients do not introduce an error in the solution using the FE method in this case. I must note that I am biased toward a FE solution, as this is the field I come from.

In your case, the coefficients also change suddenly in time. This may introduce a numerical error if not dealt with properly. You can either reduce the time step drastically after each coefficient change (this may be costly if there are a lot of such changes), or perhaps use space-time finite elements - I am not sure but I think in the space-tame case the FE scheme would not be sensitive to time jumps in the coefficients as well.

There are also other approaches to dealing with the accuracy problem due to sudden changes in the domain.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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