Hi, I'm solving a non linear diffusion equation : div(c(u) grad(u)) = 0 I'm doing that quite naively by iterating the resolution of: c(u) Laplacian(u) + grad(c(u)).grad(u) = 0 holding c(u) constant at each iteration (basically, the solution of the previous iteration) and solving the linear system via finite differences.
I saw that overall the residual decreases at each iteration, but it does so non-monotonically which makes me wonder about the convergence of such a scheme. Is this scheme proven to converge ? (is there any reference if so?). If not, is there a better simple scheme, ideally using finite differences ?
> Hi, > I'm solving a non linear diffusion equation : > div(c(u) grad(u)) = 0 > I'm doing that quite naively by iterating the > resolution of: > c(u) Laplacian(u) + grad(c(u)).grad(u) = 0 > holding c(u) constant at each iteration (basically, > the solution of the > previous iteration) and solving the linear system via > finite differences.
> I saw that overall the residual decreases at each > iteration, but it does > so non-monotonically which makes me wonder about the > convergence of such > a scheme. Is this scheme proven to converge ? (is > there any reference if > so?). If not, is there a better simple scheme, > ideally using finite > differences ?
>> Hi, >> I'm solving a non linear diffusion equation : >> div(c(u) grad(u)) = 0 >> I'm doing that quite naively by iterating the >> resolution of: >> c(u) Laplacian(u) + grad(c(u)).grad(u) = 0 >> holding c(u) constant at each iteration (basically, >> the solution of the >> previous iteration) and solving the linear system via >> finite differences.
> In cartesian coordinates and one dimension, the > expression above is usually discretized via
Thank you very much, I implemented this same discretization just before as well. Actually, the only difference I have is that I computed the midpoints as (c(u)_(i+1)+c(u)_i)/2 instead of (c(u_(i+1)))+c(u_i))/2. I don't know if this changes much.
However, I am mainly interested in knowing if the method of holding c(u) based on the previous iteration and solving as a linear diffusion equation, and iterating, is supposed to converge or not. ie. solving d/dx(c(u^{N-1})*du^{N}/dx) iterating over N, and computing the solution u^{N} at the Nth iteration.
With the above discretization (with the midpoints), on an RGB image, it converges properly on the green and blue channels, but starts oscillating around a quite moderate residual in the red channel, for some c(u). For others c(u), everything oscillate and the error remains very high...
>Hi, >I'm solving a non linear diffusion equation : >div(c(u) grad(u)) = 0 >I'm doing that quite naively by iterating the resolution of: >c(u) Laplacian(u) + grad(c(u)).grad(u) = 0 >holding c(u) constant at each iteration (basically, the solution of the >previous iteration) and solving the linear system via finite differences.
>I saw that overall the residual decreases at each iteration, but it does >so non-monotonically which makes me wonder about the convergence of such >a scheme. Is this scheme proven to converge ? (is there any reference if >so?). If not, is there a better simple scheme, ideally using finite >differences ?
I assume from your later remarks that dimension is 2. for simplicity I assume c(u) is scalar (otherwise anything is simlar, but a bit more complicated to write it down here) Hopefully c(u) >= c0>0 for all u You didn't exactly say how you chose the iteration indices in the first order product. (Thorsten is right in proposing the symmetric scheme for discretizing this selfadjoined form, you lost accuracy with your onesided meanvalue approximation) But for your current question this is irrelevant. I did what you did (I hope i am right) and got the following pde: c scalr, as said
c(u)( u_{xx} + u_{yy} ) + c'(u)( u_x^2+u_y^2) = 0 (plus your boundary conditions) now we introduce the iterated functions u^{K} : c(u^{k})*( u^{k+1}_{xx}+ u^{k+1}_{yy}) + c'(u^{k})*( u^{k}_x*u^{k+1}_x + u^{k}_y+u^{k+1}_y ) =0
and discretize this on the grid in the standard form. This finally gives something like
G(vec_u^{k}) vec_u^{k+1} = B B is the vector obtained from the boundary conditions. (*)
vec_u is the vector of the u-values on the grid. G is a matrix and will have the following form :
diag(c(u^{k}_{num(i,j))}) A + diag(c'(u^{k}_{num(i,j)})*diag(u^{k}_{num{i+1,j}-u^{k}_num(i-1,j)})*B_x + diag(c'(u^{k}_{num(i,j)})*diag(u^{k}_{num{i,j+1}-u^{k}_num(i,j-1)})*B_y
where A represents the standard five point discretization and B_x is the discretization operator for (d/dx) and similar for B_y.
Now for (*) the Banach fixed point theorem must be applied:
vec_u^{k+1} = inverse(G(u^{k})*B = Phi(u^{k})
and the next step (not that hard as it seems) is to compute the Jacobian of Phi. differentiation is to be done row and columnwise, hence essentially a scalar one. It holds for a scalar s
inserting this you see that the current u^{k} enters also the Jacobian and it remains to discuss this one.
Badly enough you see that even c''(u) will enter this calculation. But the dominant terms will depend on c' , I guess. It remains to discuss the norm of the Jacobian, and obviously this will not be less than one always.
> I assume from your later remarks that dimension is 2.
Hi, Thank you very much for this reply! Yep, the dimension is 2.
> for simplicity I assume c(u) is scalar (otherwise anything is simlar, but a bit > more complicated to write it down here)
and scalar as well :) Actually, it comes from the minimization of integral(phi(||grad(u)||)) giving c(u) = phi'(||grad(u)||)/||grad(u)||.
> Hopefully c(u) >= c0>0 for all u > You didn't exactly say how you chose the iteration indices in the first order > product. > (Thorsten is right in proposing the symmetric scheme for discretizing this > selfadjoined form, you lost accuracy with your onesided meanvalue approximation)
I now indeed use Thorsten's symmetric discretization.
> This finally gives something like
> G(vec_u^{k}) vec_u^{k+1} = B B is the vector obtained from the boundary conditions. > (*)
> Now for (*) the Banach fixed point theorem must be applied:
> vec_u^{k+1} = inverse(G(u^{k})*B = Phi(u^{k})
I don't recall Banach fixed point theorem. Wikipedia seems to say that it just states that any contraction has a fixed point, is-it just that?
> and the next step (not that hard as it seems) is to compute the Jacobian > of Phi. differentiation is to be done row and columnwise, hence essentially a scalar > one. It holds for a scalar s
> inserting this you see that the current u^{k} enters also the Jacobian > and it remains to discuss this one.
> Badly enough you see that even c''(u) will enter this calculation. > But the dominant terms will depend on c' , I guess. > It remains to discuss the norm of the Jacobian, and obviously this will > not be less than one always.
Do you know of a better scheme which always converges for any c(u) ?
Matlab's help on nonlinear FEM gives a scheme for a FEM version of the nonlinear diffusion problem. However I don't understand it from the start since it doesn't seem to be any problem for them to assemble a nonlinear pde problem (!). I just don't get how they obtain their stiffness matrix K independently of U, leading to a purely linear system in the computation of rho(u) :s And at the end, if I follow well, if I just use their simple version U^(n+1) = inv(K).F, it seems that it is doing exactly the same thing as my version using finite differences to compute K... isn't it ? Then I guess it should converge ? I get confused....
>Peter Spellucci wrote: > > I assume from your later remarks that dimension is 2.
>Hi, >Thank you very much for this reply! >Yep, the dimension is 2.
>> for simplicity I assume c(u) is scalar (otherwise anything is simlar, but a bit >> more complicated to write it down here)
>and scalar as well :) >Actually, it comes from the minimization of integral(phi(||grad(u)||)) >giving c(u) = phi'(||grad(u)||)/||grad(u)||.
>> Hopefully c(u) >= c0>0 for all u >> You didn't exactly say how you chose the iteration indices in the first order >> product. >> (Thorsten is right in proposing the symmetric scheme for discretizing this >> selfadjoined form, you lost accuracy with your onesided meanvalue approximation)
>I now indeed use Thorsten's symmetric discretization.
>> This finally gives something like
>> G(vec_u^{k}) vec_u^{k+1} = B B is the vector obtained from the boundary conditions. >> (*)
>> Now for (*) the Banach fixed point theorem must be applied:
>> vec_u^{k+1} = inverse(G(u^{k})*B = Phi(u^{k})
>I don't recall Banach fixed point theorem. Wikipedia seems to say that >it just states that any contraction has a fixed point, is-it just that?
>> and the next step (not that hard as it seems) is to compute the Jacobian >> of Phi. differentiation is to be done row and columnwise, hence essentially a scalar >> one. It holds for a scalar s
>> inserting this you see that the current u^{k} enters also the Jacobian >> and it remains to discuss this one.
>> Badly enough you see that even c''(u) will enter this calculation. >> But the dominant terms will depend on c' , I guess. >> It remains to discuss the norm of the Jacobian, and obviously this will >> not be less than one always.
>Do you know of a better scheme which always converges for any c(u) ?
>Matlab's help on nonlinear FEM gives a scheme for a FEM version of the >nonlinear diffusion problem. However I don't understand it from the >start since it doesn't seem to be any problem for them to assemble a >nonlinear pde problem (!). I just don't get how they obtain their >stiffness matrix K independently of U, leading to a purely linear system >in the computation of rho(u) :s >And at the end, if I follow well, if I just use their simple version >U^(n+1) = inv(K).F, it seems that it is doing exactly the same thing as >my version using finite differences to compute K... isn't it ? Then I >guess it should converge ? I get confused....
<< banach fixed point theorem: (in the variant you will need it, this will be no global contraction, for sure!): let D0 be convex and compact. assume Phi is C^1 on the interior of D0. If there is some norm for which the Jacobian J_Phi satisfies
||J_Phi(x)|| <= L < 1 for every x in the interior of D0.
and some y0 in the interior of D0 such that the ball
B= {z: ||z-y0||<= ||Phi(y0)-y0||/(1-L) } is contained in D0
then Phi has an unique fixed point in D0, which is also in B, and for every x0 in B the iteration x(k+1)=Phi(x(k)) converges to this fixed point at least linearly with a rate of L (everything measured in the norm ||.|| chosen for the Jacobian).
I don't have matlab/fem available, but I guess they use something like the "simplified newtons method" (the one which evaluates the Jacobian only once) in order to solve the nonliner system, but clearly they must treat the nonlinear system as such. But this is also not always convergent: you need a "sufficiently good initial guess" (and the Jacobian sufficiently accurate) for this.
Without strong conditions on c(.) there are no globally convergent schemes
But if c(.) is strictly positive (or negative) everywhere you could go back to your integral, discretize this one directly (then you get a uniformly convex functional on R^n) and then apply a globally convergent minimizer (for example cg with preconditioner, for example the Cholesky decomposition of your first stiffness matrix) which are available for this case (this includes then a stepsize scheme for possible damping the correction etc)
In this case the true full Newton's method _with damping strategy_ would also be globally convergent, but this is clearly much more ambitious than the scheme you are using presently (and which is not uncommon in this situation)
> I don't have matlab/fem available, but I guess they use something like the > "simplified newtons method" (the one which evaluates the Jacobian only once) > in order to solve the nonliner system, but clearly they must treat the nonlinear > system as such. > But this is also not always convergent: you need a "sufficiently good initial > guess" (and the Jacobian sufficiently accurate) for this.
> Without strong conditions on c(.) there are no globally convergent schemes
> But if c(.) is strictly positive (or negative) everywhere you could go back to > your integral, discretize this one directly (then you get a uniformly convex > functional on R^n) and then apply a globally convergent minimizer > (for example cg with preconditioner, for example the Cholesky decomposition > of your first stiffness matrix) > which are available for this case (this includes then a stepsize > scheme for possible damping the correction etc)
> In this case the true full Newton's method _with damping strategy_ > would also be globally convergent, but this is clearly much more ambitious > than the scheme you are using presently (and which is not uncommon in this > situation)
yeah thanks! I have no background in nonlinear pde (at school everything ended up being linear pde with a change of variable!).. but I know a bit of nonlinear optimization... Is-it a traditional way of treating nonlinear pde as nonlinear optimization ? (I mean, not directly on the integral to be minimized as you say, but on the residual of the pde).
Nicolas Bonneel <nbonn...@cs.ubc.ca> writes: > and scalar as well :) > Actually, it comes from the minimization of integral(phi(||grad(u)||)) > giving c(u) = phi'(||grad(u)||)/||grad(u)||.
>> and scalar as well :) >> Actually, it comes from the minimization of integral(phi(||grad(u)||)) >> giving c(u) = phi'(||grad(u)||)/||grad(u)||.
> What is phi?
> Nicolas
Hi, I was planning to experiment with several phi... * phi(||grad u|| = sqrt(1+grad(||u||)^2) * phi(||grad u|| = ||u||^n (n=2 corresponds to the linear pde) * gaussian-like * any exponential-like function * maybe others.... :p
Nicolas Bonneel <nbonn...@cs.ubc.ca> writes: >> Hi, >> I was planning to experiment with several phi... >> * phi(||grad u|| = sqrt(1+grad(||u||)^2) > phi(||grad u|| = sqrt(1+||grad u||^2)
This is the functional for the (minimal surface) Plateaux problem.
>> * phi(||grad u|| = ||u||^n (n=2 corresponds to the linear pde) > phi(||grad u|| = ||grad u||^n
And this is the functional for the p-Laplacian.
There is already a lot in the literature about approximation of the corresponding PDEs, especially in connection with finite elements. You can also assume that (most of) the finite element results have their counterpart for structured meshes and finite difference discretizations.
|
