How to solve filtration equation in porous media

Hi every body

I am considering to solve a PDE which has the form:

dC/dt = d2C/dx2 + dC/dx + C    x is from 0 to L

Boundary condition:

dC/dx = 0 @ x = L;
C(0,0)= Co;
C(0,x)= 0;
C(t,0)= 0

Can Matlab solve such a equation like that? If it can, how can I set it up?

Thanks for any help.

Answer 1

[Grzegorz Knor]

pdepe funtion can solve it. Simple example:

http://www.mathworks.com/matlabcentral/answers/16550-pdepe-function

[NGOC HONG]

I looked at that already but my equation is quite different. Mine has two additional terms dC/dx and C. I also notice that Matlab also has a function that can solve equations in the form dC/dt = d2C/dx2 + dC/dx but I am still not sure I will work for mince case.

Thanks

Answer 2

[Walter Roberson]

No, MATLAB cannot solve that. The boundary conditions you give are unsolvable unless Co happens to be 0.

You have set C(0,x) = 0 and C(t,0) = 0, so C(0,0) has to satisfy both those conditions which requires that C(0,0) = 0 which contradicts your requirement that C(0,0) = Co unless Co happens to be 0.

You might wish to refer to an older thread in which I show how to work with some boundary conditions with dsolve(): http://www.mathworks.com/matlabcentral/answers/18495-ordinary-differential-equation-with-boundary-condition-in-z-0-and-in-z-l

The bad news is that your conditions are only solvable if C(t,x) = 0 as a constant. Your requirement that C is 0 along both axes is only consistent with your original function when C is 0 everywhere.

If you remove both requirements that C(0,x) = 0 and C(t,0) = 0 but take in to account dC/dx = 0 @ x = L then Maple's pdsolve() does offer an additional answer to

pdsolve([diff(C(t, x), t) = diff(C(t, x), x, x)+diff(C(t, x), x)+C(t, x), eval(diff(C(t, x), x), x = L) = 0], C(t, x));

namely,

C(t,x) = _C1 * _C2 * exp(_c[1]*t-(1/2)*x) * ((-1+(-3+4*_c[1])^(1/2)) * exp((1/2)*(-3+4*_c[1])^(1/2)*(2*L-x)) + exp((1/2)*x*(-3+4*_c[1])^(1/2)) * (1+(-3+4*_c[1])^(1/2))) / (1+(-3+4*_c[1])^(1/2))

where _C1 and _C2 I am sure are constants of integration, but I am not sure what _c[1] is intended to represent (another constant of integration would be most likely.)

When I work this possibility through to try to add in either the boundary condition C(0,x)=0 or C(t,0)=0, the solutions again reduce to C(t,x)=0 always. However, if those were mistaken boundary conditions, then the above gives a possibility for other solutions.

Unfortunately MuPAD does not have a pdsolve() function, but possibly you could work with a boundary-value solver such as bvp4c() after manually converting your pde to a set of ode. Might be tricky because one of your boundary positions is symbolic but the bvp*c series are finite difference solvers.

[NGOC HONG]

If I change the boundary condition into: C(0,x)=0, C(0,0)=0, C(1,0)=Co, and dC/dx = 0 @ x = L then it will be solvable with Matlab, right?

[Walter Roberson]

If you examine the C(t,x) expression above, and substitute 0 for t in order to examine the boundary C(0,x), you end up with the multiplication of 3 terms. If you then study the conditions under which the resulting expression can be 0, you end up with three possibilities: either _C1 is 0, or _C2 is 0, or a particular expression in _c[1] is 0 (the expression has zeros at _c[1] = 3/4 exactly, and at _c[1] = 0.98 approximately.) But no matter which one of those three terms you set to 0 as a constant, you end up driving C(t,x) to 0 always for all t and x.

You can work it the other way, too, letting x be 0 in the C(t,x) expression shown above; you come out with the same conditions for the expression to be 0.

Your choice of the symbolic value Co as a boundary value prevents you from using bvp4c or bvp5c as your boundary value solver, as those can only work with boundaries that are purely numeric. There does not appear to be a MuPAD PDE solver, but if you were to convert your PDE to a set of ODE then you could use the MuPAD dsolve() to come out with the basic equations, after which you would have to do some work in order to solve for the boundary conditions. I do not know at the moment if it is feasible to automate that boundary condition solving in MuPAD. I guess it must be possible or else Maple would not be able to do it, but I do not have a feeling for how much work it would involve. (One would tend to think that if it was easy, MuPAD would already have an implementation, for ease of use if nothing else.)

[Walter Roberson]

If there is a time tn for which C(tn,0) = 0, then the only solution is C(t,x) = 0 always. This includes the condition C(0,0)=0 as a subset. Likewise if there is an x, xm, for which C(0,xm) = 0, then the only solution is C(t,x) = 0 always.

Thus, imposing any boundary condition of 0, along either the t or x axis, forces the entire solution to the constant 0.

Imposing a non-zero boundary condition along either the t or x axis does NOT force the entire solution to be the constant 0.

Imposing an non-zero boundary condition along the t axis and another non-zero boundary condition along the x axis, appears to be self-consistent, but the expression for it gets very long:

If the boundary condition at C(m,0) = Cm0, and the boundary condition at C(0,n) = C0n, then the constant _c[1] shown in the above Answer, is forced to be,

(1/4) / L^2 *(3*L^2 + RootOf(L^2*C0m^2*exp(_Z*(2*L+m)/L) - 2*L^2*C0m*Cn0*exp((1/4)*(8*_Z*L^2+2*_Z*L*m+3*t*L^2+_Z^2*t-2*m*L^2)/L^2) + L^2*Cn0^2*exp((1/2)*(4*_Z*L^2+3*t*L^2+_Z^2*t-2*m*L^2)/L^2) - _Z^2*C0m^2*exp(_Z*(2*L+m)/L) + 2*_Z^2*C0m*Cn0*exp((1/4)*(8*_Z*L^2+2*_Z*L*m+3*t*L^2+_Z^2*t-2*m*L^2)/L^2) - _Z^2*Cn0^2*exp((1/2)*(4*_Z*L^2+3*t*L^2+_Z^2*t-2*m*L^2)/L^2) - 2*L^2*C0m^2*exp(_Z*(L+m)/L) + 2*L^2*C0m*Cn0*exp((1/4)*(4*_Z*L^2+6*_Z*L*m+3*t*L^2+_Z^2*t-2*m*L^2)/L^2) + 2*L^2*C0m*Cn0*exp((1/4)*(4*_Z*L^2+2*_Z*L*m+3*t*L^2+_Z^2*t-2*m*L^2)/L^2) - 2*L^2*Cn0^2*exp((1/2)*(2*_Z*L^2+2*_Z*L*m+3*t*L^2+_Z^2*t-2*m*L^2)/L^2) - 2*_Z^2*C0m^2*exp(_Z*(L+m)/L) + 2*_Z^2*C0m*Cn0*exp((1/4)*(4*_Z*L^2+6*_Z*L*m+3*t*L^2+_Z^2*t-2*m*L^2)/L^2) + 2*_Z^2*C0m*Cn0*exp((1/4)*(4*_Z*L^2+2*_Z*L*m+3*t*L^2+_Z^2*t-2*m*L^2)/L^2) - 2*_Z^2*Cn0^2*exp((1/2)*(2*_Z*L^2+2*_Z*L*m+3*t*L^2+_Z^2*t-2*m*L^2)/L^2) + L^2*C0m^2*exp(_Z*m/L)-2*L^2*C0m*Cn0*exp(-(1/4)*(-6*_Z*L*m-3*t*L^2-_Z^2*t+2*m*L^2)/L^2) + L^2*Cn0^2*exp(-(1/2)*(-4*_Z*L*m-3*t*L^2-_Z^2*t+2*m*L^2)/L^2) - _Z^2*C0m^2*exp(_Z*m/L) + 2*_Z^2*C0m*Cn0*exp(-(1/4)*(-6*_Z*L*m-3*t*L^2-_Z^2*t+2*m*L^2)/L^2) - _Z^2*Cn0^2*exp(-(1/2)*(-4*_Z*L*m-3*t*L^2-_Z^2*t+2*m*L^2)/L^2)) ^2)

In the above, the long expression RootOf involving the variable _Z indicates that the value of the RootOf expression is the set of values _Z for which the expression inside the RootOf() becomes zero -- i.e., the value returned but the construct is the roots with respect to _Z of the expression inside. If you could manage to get a nice formatting of the expression, you would see that the expression is a mix of terms: constant expressions in the boundary positions (L, m, and n), and the boundary values (C0m, and Cn0), together with terms in which _Z^2 appears as a multiplier, and exp() terms in which _Z to the first power appears in combination with L, m, and t.

Unfortunately, even if you were to substitute particular numeric values for the boundary conditions, the symbolic solution would be a mess, and would probably be quite uncommon that it could be solved down to explicit formula: you would probably usually need to find numeric roots in _Z.

[NGOC HONG]

Sorry for making it more complicated. I checked the boundary conditions and the correct one was:

C=Co @ x=0,0=<t=<sigma(t). sigma(t) means very small t.

C=0 @ x=0, t>sigma(t).

dC/dx=0 @ x=L

Is is solvable now?

Thanks for you patient!

[Walter Roberson]

In order for C(0,t) to change abruptly from C0 to 0, C(0,t) would have to involve a close cousin of the Heaviside step function, differing from MATLAB's symbolic Heaviside in the value at exactly sigma(t). MATLAB (and other sources) usually define Heaviside(sigma(t)) as being 0 at x < sigma(t), and 1 at x > sigma(t), and 1/2 at sigma(t) exactly. The 1/2 at sigma(t) is a problem. Some sources define the the value at the exact parameter to be undefined (or user defined) instead of 1/2. You could use MATLAB's symbolic piecewise() to implement a 0/1/1 scheme, but then you would have trouble differentiating it.

You could remove this problem if you could allow C = C0/2 at x=sigma(t) .

I would have to examine further to see if it would be possible using Heaviside.