Timeline for Numerical simulation of a reaction-diffusion system on MATLAB with finite difference discretization of spatial derivative
Current License: CC BY-SA 3.0
16 events
| when toggle format | what | by | license | comment | |
|---|---|---|---|---|---|
| Mar 3, 2016 at 6:00 | comment | added | Chris Rackauckas | Oh yeah, if you're using an adpative solver don't worry about changing dt to check the order of accuracy. You can just assume your temporal integrator (ode15s) is correct. I forgot about that and so for example if you coded your own RK4, you would want to check you get 4th order accuracy. | |
| Mar 2, 2016 at 18:50 | vote | accept | WYSIWYG | ||
| Mar 2, 2016 at 18:50 | comment | added | WYSIWYG |
Thanks. I was trying it on a Turing problem. I'll do the basic problems too. Diffusion seems to work as I can see the concentration becoming uniform. I have to check the dynamics though. Reaction works like well mixed system (ODE) when diffusion is set to zero. I didn't get your point on halving and doubling dt. I do not set dt. I use ode15s to solve the spatially discretized problem. Do you mean tolerance? I can change dx and see if the solution is similar.
|
|
| Mar 2, 2016 at 15:02 | comment | added | Chris Rackauckas | Something like this or this. If you have the right linear behavior and can reproduce non-trivial nonlinear behavior that's a good sign it's working. | |
| Mar 2, 2016 at 15:00 | comment | added | Chris Rackauckas | First try the diffusion equation (no reaction). There is a known solution via Fourier transforms that you can test against. You should check that your order of accuracy is 2 (evaluate by halving/doubling dx a few times and graph it). Then set diffusion to zero and test a reaction equation. Usually I would use some linear reaction equation since there is an actual solution, and you should get the order that matches your method (halving and doubling dt this time). When you have all of that in order, I usually make sure it can solve a Turing patterning problem. | |
| Mar 2, 2016 at 7:30 | comment | added | WYSIWYG | I have done some modifications in the code such that there are no loops and the computations are done on the vectorized/matrix forms. I also split the diffusion from the reaction terms as suggested in Randy LeVeque's book. The speed has tremendously increased. However, I need to know if this program is actually solving the problem correctly. Can you suggest me some simple published solutions that I can replicate? Do you think it is a good idea to share my program somewhere for a review (Code Review??) | |
| Feb 12, 2016 at 14:58 | comment | added | Chris Rackauckas | This vectorized form in MATLAB should be by far the fastest. However, there are some details in the implementation as to why it's not as fast as possible in more efficient langauges I have a blog post about it if you're interested, so if the vectorized version is too slow, then take look at writing your function as a MEX (C) call or in Julia in the devectorized form. Or you can go all the way out to Fortran, but that's a little inconvenient for graphing. | |
| Feb 12, 2016 at 14:58 | comment | added | Chris Rackauckas | However, given the way de-referencing works, if you have lots of reactions, it could be more efficient to just save the matrix of u(1:26:end),u(2:26,end), etc., and the act of building this matrix in its most efficient form is reshape. | |
| Feb 12, 2016 at 14:51 | comment | added | Chris Rackauckas | I would be highly surprised if that's slower than a loop. Another thing you can do is in your vectorized code you can call "sub-vectors of the reactants". For example, if your vector is all the reactants at space j, then all the reactants at j+1, etc., then you can get the vector of all of the reactant 1 via u(1:26:end). So you can add at all points in space for two reactants via u(1:26:end)+u(2:26:end). If you grouped by space, you would just make the middle number the number of space points. This is pretty much what your loop is likely doing. | |
| Feb 12, 2016 at 10:18 | comment | added | WYSIWYG | Yes I can reshape the vector to a matrix and apply f on all the columns and again reshape back to a vector. This, I guess, would take more time than running a loop. | |
| Feb 12, 2016 at 10:13 | comment | added | WYSIWYG | No I didn't mean time. We can perhaps discuss this better in chat if you are free for that. I just meant that f takes all reactants at a given space and time point. If u is a long vector (m×n,1) of all reactants(n) at all spatial points(m) either listed reactants first or timepoints first. The f should be redefined such that it operates on this huge vector; basically reapply the same reaction rules for all the m-spatial points. This will require a loop. | |
| Feb 12, 2016 at 4:45 | comment | added | Chris Rackauckas | $f$ is applied to all reactants at all points in time? That makes no sense in just about any model. I edited for clarity that $u_i$ is the vector of all reactants at the time point $i \Delta t$. In MATLAB you write $f$ as a vector valued function that takes in all of the reactants at time $i$ and spits out $f$ evaluated at each point in space (you may need to reshape the vector inside of $f$ and use a bunch of .* and ./). Then the code is simply has u(i+1)=A*u(i) + f(u(i)) as the update function with no loop. Let me know if you need more details. | |
| Feb 12, 2016 at 4:40 | history | edited | Chris Rackauckas | CC BY-SA 3.0 |
added 51 characters in body
|
| Feb 12, 2016 at 4:32 | comment | added | WYSIWYG | I don't need a stiff solver because the rates are more or less similar for all variables. Like you said I can represent the system as you showed above. However, the equation is more like this $u_i''(x,t)=Au_i+f(u)$ i.e. the nonlinear reaction function requires all $u_i$. So it is difficult to use the matrix form. I'm thinking of better ways which I think I'll optimise as I keep working. My current issue is to understand if the approach is good enough. | |
| Feb 12, 2016 at 2:45 | history | edited | Chris Rackauckas | CC BY-SA 3.0 |
edited body
|
| Feb 12, 2016 at 2:38 | history | answered | Chris Rackauckas | CC BY-SA 3.0 |