[R] finite difference scheme for 2D differential equations

[Thomas Petzoldt]

Hello Subodh Acharya,

I've been away for a field trip for two weeks, and I guess you have 
already found a solution for your problem.

If your question is still open, you may consider to look into package 
deSolve that now provides functions (ode1D, ode2D and ode3D) for solving 
1-3 dimensional equations using the method of lines approach. The 
package comes with many examples and extensive documentation and there 
is also a JSS publication:

http://www.jstatsoft.org/v33/i09

If you need more help, please write to the dedicated mailing list for 
dynamic models:

https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

Hope it helps

Thomas Petzoldt

PS: You may also have a look into package ReacTran.


At 26.04.2010 03:09 Subodh Acharya wrote:
> Hello everyone,
>     I am trying to solve 2D differential equations using finite difference
> scheme in R. I have been able to work with the equations with only one
> spatial dimensions but I want to extend it to the two dimensional problem.
> For example i can simulate one dimensional diffusion using a code like the
> following. But I want to write a similar code for,say, a two dimensional
> diffusion equation. Any kind of help/advice is highly appreciated.
>
> Here is how I did for the 1D equation (in this case for 1D diffusion)
> Equation: dc/dt  = D*d^2c/dx^2
>
>
> dx = 10
> dt = 0.1
> D = 15
> n = 20
> tstep = 200
> C = 50
> dt = 0.1
> dx = 10
> conc<-matrix(0,tstep,n)
> conc[1,1:10] = C
>   for (i in 2:tstep) {
> for (j in 1:n){
> if (j==1){
> conc[i,j] = conc[i-1,j] - dt*(D/dx^2)*(conc[i-1,j] - conc[i-1,j+1])
> }
> conc[i,j] = conc[i,j]
>   if(j>1&&  j<  n){
> conc[i,j] = conc[i-1,j] + dt*(D/dx^2)* (conc[i-1,(j-1)] - 2*conc[i-1,j] +
> conc[i-1,(j+1)])
> }
> if (j==n){
> conc[i,n] = conc[i-1,n] + dt*(D/dx^2)* ( conc[i-1,n-1] - conc[i-1,n])
> }
> conc[i,j] = conc[i,j]
> }
> }
>
> Now in 2D the equation will be like this
>
> dc/dt  = Dx*d^2c/dx^2 + Dy*d^C/dy^2
>
> So that when I solve it I will get C(x, y, t) at each node of the grid.
>
> Thanks for the help in advance.