fast fourier transform of classical diffusion equation

[kedaR]

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

[nirgrahamuk]

@kedaR

{
c[i] =

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.

[kedaR]

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

[nirgrahamuk]

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]] <- 1

This will succeed , but abc is a list

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

This will also succeed , but abc is a numeric vector

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

[kedaR]

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

[nirgrahamuk]

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)

[kedaR]

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