# fast fourier transform of classical diffusion equation

My `version` output in RStudio,

``````platform       x86_64-w64-mingw32
arch           x86_64
os             mingw32
system         x86_64, mingw32
status
major          4
minor          0.4
year           2021
month          02
day            15
svn rev        80002
language       R
version.string R version 4.0.4 (2021-02-15)
nickname       Lost Library Book
``````

My MRE:

``````N = 32;
dx = 1.0;
m = 1;

for (i in 1:5) {
c[i] = 0.5*(1+sin(2*pi*m*i*dx/N));
}

plot(c)
par(new=TRUE)
halfN = N/2;
delk = 2*pi/N;
dt = 0.5;

for (m in 1:80) {
chat = fft(c);
for (i in 1:N) {
if ((i-1) <= halfN) k <-
(i-1)*delk;
if ((i-1) > halfN) k <-
(i-1-N)*delk;
k2 = k*k;
chat[i] = chat[i]/(1 + k2*dt);
}
c = fft(fft(c), inverse = TRUE)/length(c);
}

plot(c)
``````

I get an empty output plot as follows: I am expecting to get the sinusoidal wave as initialized in the 1st for loop & `fft()` of the same to solve the diffusion equation.

Where am I going wrong with the code

what is this ?
`c` is a name of a very useful function that you should avoid overloading, and you need to define variables that you will access the indexes of.

1 Like

Oh, yes. That is so silly of me!

``````N = 32;
dx = 1.0;
m = 1;
for (i in 1:5) {
C(i) = 0.5*(1+sin(2*pi*m*i*dx/N));
}
plot(C);
par(new=TRUE)
halfN = N/2;
delk = 2*pi/N;
dt = 0.5;

for (m in 1:80) {
Chat = fft(C);
for (i in 1:N) {
if ((i-1) <= halfN) k <-
(i-1)*delk
if ((i-1) > halfN) k <-
(i-1-N)*delk
k2 = k*k;
Chat(i) = Chat(i)/(1 + k2*dt);
}
C = fft(fft(C), inverse = TRUE)/length(C);
}
plot(C)
``````

Thank you @nirgrahamuk. I changed the variable to `C` now.

In the following chunk,

``````for (i in 1:5) {
C(i) = 0.5*(1+sin(2*pi*m*i*dx/N));
}
``````

I am still receiving the error,

``````Error in C(i) <- 0.5 * (1 + sin(2 * pi * m * i * dx/N)) :
object of type 'closure' is not subsettable
``````

ok, theres same issue again as stats::C is a function also...
but aside from that, as I said in my previous post, you cant expect to get at the indexes of objects you havent created yet. consider the following example
This will fail

``````abc[] <- 1
``````

This will succeed , but abc is a list

``````abc <- list()
abc[] <- 1
``````

This will also succeed , but abc is a numeric vector

``````abc <- vector(length=5)
abc[] <- 1``````
1 Like

Hello @nirgrahamuk,

``````N = 32;
dx = 1.0;
m = 1;
for (i in 1:5) {
a[i] = 0.5*(1+sin(2*pi*m*i*dx/N));
}
plot(a);
par(new=TRUE)
halfN = N/2;
delk = 2*pi/N;
dt = 0.5;

for (m in 1:80) {
ahat = fft(a);
for (i in 1:N) {
if ((i-1) <= halfN) k <-
(i-1)*delk
if ((i-1) > halfN) k <-
(i-1-N)*delk
k2 = k*k;
ahat[i] = ahat[i]/(1 + k2*dt);
}
a = fft(fft(a), inverse = TRUE)/length(a);
}
plot(a)
``````

This MRE compiles without any error and gives the output a plot,

Thank you

I think you may be using automatic save and restore of your workspace. this is not recommended for reproducible work.

again, in the section

``````for (i in 1:5) {
a[i] = 0.5*(1+sin(2*pi*m*i*dx/N));
}
``````

a[i] is accessed, but a has not been defined as anything appropriate and would be expected to get the following error

``````Error: object 'a' not found
``````

Therefore the code you shared will not work, unless by accident the user (including yourself) has an appropriate a object that can be indexed etc.

There is an easy fix, it is to define a in advance of accessing its parts a[i] i already gave examples so I wont repeat here.

Glad you are happy with your script..., but you are in danger of being 'caught out' if down the road when your r session is in another state (lacking an `a` for example), at which time your previously working code will fail)

1 Like

My fruitful discussions with both of you. Asked the same Q on StackOverflow as well here Initialized the vector `a` as

``````a = vector("numeric", 5)
``````

before calling/declaring it in the `for` loop.

Thank you for your valuable inputs and for letting me in on my errors This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.