[R] finite difference scheme for 2D differential equations
Thomas Petzoldt
thpe at simecol.de
Wed May 12 08:25:57 CEST 2010
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.
More information about the R-help
mailing list