Solving a stiff PDE system with operator splitting

147 Views Asked by At

I have a particular system of PDE that have given me a lot of grief and I wanted to try to solve it using operator splitting however I don't fully understand how its done. I'll start by explaining what I've done and my issues with it

The equations of interest are

(1) $ u_t=iu_{xx}+f(\epsilon)u$

(2) $\epsilon_t=q(x)-(1+g(u))\epsilon$

(3) $u(0,t)=u(L,t)=0, u(x,0)=p(x), \epsilon(x,0)=0$

Where the function q(x) is positive over the domain.

My basic attempt is to apply finite difference on the derivative and solve the system of equations directly. I used backwards Euler and treated the nonlinear part of (1) explicitly. Equation (2) can be solved exactly in backwards Euler (and in crank-nicholson as well).

(4) $\epsilon^{n+1}=\frac{\epsilon^n+\Delta tq(x)}{1+\Delta t(1+g(u^{n+1}))}$

So my general plan, was to solve equation (1), feed those results into equation (4), and repeat. While it does work it has a major flaw in that when I change the time step size the resultant solutions distort.

I want to move on to operator splitting but I don't understand how to apply it.

Splitting equation (1) is simple

(5) $u_t=(A+B)u$

(6) $A=\frac{\partial^2}{\partial x^2}, B=f(\epsilon)$

Equation (2) I'm not sure how to split since it's inhomogeneous. Would it go like this?

(7) $\epsilon_t=q(x)-(C+D)\epsilon$

(8) $C=1, D=g(u)$

Or could I say the operator C is

(9) $C\epsilon=q(x)-\epsilon$

Theoretically the solutions for operator splitting would be

(10) $u(t+h)=e^{(A+B)h}u(t)$

(11) $\epsilon(t+h)=e^{(C+D)h}\epsilon(t)$

If I did strang splitting I would have for instance

(12) $u(t+h)=e^{0.5hA}e^{hB}e^{0.5hA}u(t)$

But how do I actually use this? I know you can expand the exponential operators

(13) $e^{hA}\approx1+hA+\frac{(hA)^2}{2}+...$

But I don't see how this makes the new problem any more solvable.

I've seen you can break the operators into sub problems which would be very useful because the equations $u_t=Au$ and $\epsilon_t=C\epsilon$ are both exactly solvable if the boundary/initial conditions are unchanged. But again I don't see how use them.

What's the procedure here?