Solving a diffusion-reaction equation using method of lines

Hi all, I am trying to solve a diffusion-reaction equation (under cylindrical co-ordinate) using method of lines (MOL). The partial differential equation, initial and boundary conditions, and the corresponding derivations can be seen in the attachment. The corresponding code has also been uploaded. I also simulated the PDE using the pde solver in Matlab, and found that my MOL code cannot obtain the right results. Moreover, what weird is that the results of MOL depend on the number of spatial grids (the variable n in my code). I have checked my code many times and did not find any errors. Can anyone help me with this? Thanks in advance.

[J. Alex Lee]

i don't think your BC is implemented correctly...wouldn't you just want [1,0,0...] in your top row of A so that the equation just forces C1=BC? And same for last row?

[Dadi Bi]

Hi, thank you very much for having a look at my question. Sorry, I did understand your reply. What do you mean that I want the top row of A to be [1,0,0...]?

When j=1, Eq. (5) becomes

.

So in Eq. (6) of the pdf file, the first and second elements of the first row of A are and , respectively. As we can obtain due to the BC at the left boundary , the fisrt element of the first row of A will be updated to

.

Is this wrong?

[J. Alex Lee]

so you actually have n+2 nodes, and are only solving the "middle" n of them. I guess this makes more sense in your case for BC at left because you have an asserted time dependent function.

But then still confused why build your BC into the equation for the 1-th node? it should naturally be taken care of without "asserting" any BCs on your internal node FDEs, i.e., not sure you need the b vector at all?

[Dadi Bi]

Thanks! Yep, you are right. I have n+2 nodes and only solve the middle n nodes. For example, in the following figure, we have five nodes in total, and I only solve the x_1, x_2, and x_3 nodes.

As mentioned in my last reply, for the first node x_1, the concentration at this node can be expressed:

,

where the value of can be obtained through the BC at x_0. Due to the function , I think we still need the vector b.

Indeed, the derivation follows the MOL example (page 21) in this book (link: http://svmoore.pbworks.com/w/file/fetch/78638864/PDE_Numerical_Book_Olsen-Kettle.pdf). And I used the author's code as a template (has been uploaded). I feel my derivation and coding are right but do not know why I cannot obtain the same result using the pde solver.

[J. Alex Lee]

Your explanation in the comment is clear, now that I understand the BCs are flux conditions, based on your equations it is

The method of not solving for the end nodes seems clever, and your equations do look ok to me. The only comment I would make is that your flux condition is asserted using a 2 point FDE (first order accurate) whereas your other FDE's are 2nd order accurate. This might result in subtly different answers. How are you generating your "reference" solution (pdepe?) and in what way is your answer different from it?

By teh way, you can also re-institute the much better way of constructing A in the example

B = [ (s - s./(2*(1:n)')) , (-2*s-kd)*ones(n,1) , ( s + s./(2*(-1:n-2)')) ];
A = spdiags(B,-1:1,n,n);

[Dadi Bi]

Many thanks and much appreciate that you had a look at my question even on weekend. You are right that the BCs are flux conditions. It is also true that the first and second spatial derivations are second order approximations while the BCs are first order approximation. I will use the first order approximation for first order spatial derivation to see what will happen.

For the "reference" solution, I used the pdepe function (see attachment "pdex1.m"and https://uk.mathworks.com/help/matlab/math/partial-differential-equations.html). The results of PDE solver and the MOL with n=99 are shown in the following figures. As you can see, the curves using MOL has the same trend with that using PDE solver, but the concentration is different.

As I mentioned in the question, what weird is that the results of MOL are dependent on the variable "n". The larger the n, the lower the concentration.

Finally, thanks for your advice on how to construct matrix A.

[J. Alex Lee]

In your pdepe example, you are asserting your flux condition at some finite value of r rather than at r=0. If you change your boundary location to r=0, the pdepe solver seems to fail. If you change your boundary to match that of your original formulation (for larger n, you get smaller r0), you'd see a dependence on the boundary location similarly to what you see.

In hindsight, the BC at r=0 violates radial symmetry, which seems to be problematic.

[Dadi Bi]

Thank you very much! At first, I considered n=100, and thus the first node will be 0.05. I noticed that the pde solver failed when x = linspace(0,5,100) so that I just set x = linspace(0.05,5,100) in pdepe. However, by doing so, you are right the left BC became instead of .

Let’s assume that we consider and the left and right BCs are and . In such a case, I think my MOL code shoud get the same results with the pdepe example if I change x1=0 to x1=0.05 and modify dx as dx=(x2-x1)/(n+1) (dx also need to be modified in "Uprime.m"). However, their results are still very different. What do you think about this?

In the book I mentioned in the last reply, the nuclear waste example in page 43 also considers a diffusion-reaction equation under the cylindrical co-ordinate with 0<r<100, which is

.

The author solved this equation using backward Euler method (implicit) and the corresponding code can see "NuclearWaste". The results of "NuclearWaste" is independent on the number of nodes n. I also solved this equation using the pdepe solver. I set x = linspace(0,100,300) and the results of pde solver are the same with the results of backward Euler method. And I am suriprised that the pde solver works when the starting node is 0.

[J. Alex Lee]

In the documentation it is stated that for pdepe that when m=1 then any boundary condition at r=0 is ignored and the symmetry condition is assumed.

The nuclear waste has a symmetry condition at r=0, so it should be ok.

Again, when you say results are very different, are you remembering to compare only the internal nodes of the pdepe solution? Remember because you have very high gradients near r=0, that extra node is going to shoot up to a high value.

[Dadi Bi]

Thank you for pointing out the setting of BC at r=0 when m=1 in the documentation. The left BC of the nuclear waste exactly matches this setting so I understand why it works.

If we consider and the BCs and , I plot the results of pdepe, MOL with n=99 and 299.

The distance ranges for pdepe, MOL with =99, and MOL with n=299 are [0.05,5], [0.0995,4.9505], [0.0665,4.9835], respectively. Although these three curves are plotted using different nodes, I think they should have the same concentrations. Moreover, for MOL, the value of variable "n" should only be related with the concentration accuracy, and it is weird to see such a huge concentration difference. Am I right?

[J. Alex Lee]

where is your code for generating these plots? You should be comparing

sol = pdepe()
u = sol(:,:,1);

figure;
mesh (t,x(2:end-1),u(:,2:end-1)')

to within a transpose of whatever is in u.

That's assuming you can control the number of nodes in pdepe, and you should compare n=101 for pdepe to n=99 for your MOL, and n=301 for pdepe to n=299 for your MOL

[Dadi Bi]

Hi, thanks! I revised my code and uploaded them. The results are still quite distinct.

[J. Alex Lee]

Are you sure your xmesh is consistent across your solution methods? Your MOL needs to be define from x1+dr to a-dr, and the term in your denominator (2*i) needs to correspond to the elements in xmesh correctly. At first glance, i think you need to add an offset of 0.05 (or whatever x1 is) in that denominator.

Also, when setting x1 = 0.05 in a MOL, i'm seeing a dramatic difference between using first order vs second order difference equations at the boundaries. Likely the pdepe is internally using consistent order of discretization, so you'd want to use second order difference eqs at the ends if you want to compare apples to apples.

[Dadi Bi]

Hi, sorry for the late response. You are right that there should be an offset in the denominator, and I have revised my code accordingly. For the BCs, I used the first order difference. I considered . I present some results in the following. We can see a closer match between MOL and pdepe with the increase of variable n, which demonstrates the correctness of the MOL code.

For the BCs, how can I use a second order difference for them. Take the left BC for exmaple, as we calculate the concentration from point and the calculation of is dependent on , it seems that we can only use the first order difference for , that is .

By the way, when the starting node is 0.05, the results are almost the same for different n (we can observe this from the above four figures). However, I am still confused that the results are quite different for different n when the starting node is 0.

[J. Alex Lee]

the LHS of your last equation is the finite difference equation estimating the left-handed first order accurate first derivative. You can replace that with any FDE you want, in particular for a left-handed second order accurate first derivative estimate:

see, e.g., https://en.wikipedia.org/wiki/Finite_difference_coefficient.

But I still have doubts about the non-zero flux at the center-line; maybe pdepe SHOULD fail when you try to assert anything but a symmetry condition there. If you are modeling a physical system, it seems reasonable that any type of line-source/sink is anyway an idealization of a "thin embedded wire"-like source so that your strategy of nudging the left boundary a bit is actually the more sensible one.

For your mesh dependence study, what are you comparing? I would compare the value of at some specific time. If you are comparing values, are you compensating for your mesh size change by changing the location of your left boundary to ?

[Dadi Bi]

Hi, thanks for link. I revised my code accordingly and didn't see a dramatic difference between using first order and second order difference equations for the left BC. The following is some results:

For the non-zero flux , α is the production rate of a type of molecule, is the number of cells in a columnar colony with radius and depth . Although the colony is not a ponit, I consider the colony as a ponit source and study the diffusion of molecules generated by cells. When I first use pdepe and set the starting node as 0, the pdepe retured some errors (sorry I forgoe these errors), and this motivates me to set the starting node as 0.05. But now, when the starting node is 0, the pdepe works and retures a zero concentration, as shown in the following figure.

I think this result does makes sense, as pdepe automatically enforces the left BC as zero and the right BC and IC are both zero, so the concentration should be always zero. As you said that nudging the left boundary a bit is actually the more sensible one, the following question is how long is the left BC nudged. I simulated the MOL with different starting node and their results are quite different.

For the mesh dependency, I am comparing the variable "U" in MOL with different n when . Or for simplicity, I compare the maximum concentration at t=24. I noticed that when , U will rapidly converge to a non-zero concentration with the increase of n (usually, n=1099, U is quite steady). However, when ，the max concentration of U drops from 3800 when n=99 to 300 when n=3099.

[J. Alex Lee]

once you change r_0, you are solving different problems, so of course you expect different answers.

still not clear if you doing the right comparison...if you are assessing mesh dependence of MOL, then you must remember that every time you change your mesh, your value of C1 (which is largest) is NOT at r_0, but at r_0+dr, and dr is always changing.

[Dadi Bi]

Yes, for the mesh dependence, although the mesh nodes are different, I think the concentration changes for different meshes should be almost the same (I am not comparing two different nodes but comparing the whole concentration change). For exmaple, if the concentration follows a function y=2x+1, for mesh settings x_mesh1=0:0.01:2 and x_mesh2=0:0.03:2, these two different meshes can result in the same concentrations 2x_mesh1+1 and 2x_mesh2+1.

[J. Alex Lee]

If you are comparing C1 values to C1 values (the first value in your mesh), instead of C0 values to C0 values, your answers WILL depend on your mesh because of the way you are defining the mesh.