N.Kalosha at math.ru.nl
2006-Apr-20 10:26 UTC
[Rd] bug: code not working as expected (PR#8783)
Hi, I've attached two files with the sources for a function to implement the finite difference method for solving a particular differential equation. One of them works correctly, the other gives wrong results or returns an error, depending on the version of R. The difference between them is that in the 'broken' version in line 42 I check if the items in the two-dimensional array are bigger than a certain value, and in the working one I do it in a separate loop. A 2.1.1 build for Solaris returns the following error Error in if ((X - (j - 1) * dS) > f[i, j]) f[i, j] = (X - (j - 1) * dS) : missing value where TRUE/FALSE needed On Windows both the stable 2.2.1 and 2.3.0rc gui versions will silently produce incorrect data. Nikolai.
N.Kalosha at math.ru.nl writes:> This is a multi-part message in MIME format. > --------------020909040800030906040005 > Content-Type: text/plain; charset=KOI8-R; format=flowed > Content-Transfer-Encoding: 7bit > > Hi, > > I've attached two files with the sources for a function to implement the > finite difference method for solving a particular differential equation. > > One of them works correctly, the other gives wrong results or returns an > error, depending on the version of R. > > The difference between them is that in the 'broken' version in line 42 I > check if the items in the two-dimensional array are bigger than a > certain value, and in the working one I do it in a separate loop. > > A 2.1.1 build for Solaris returns the following error > > Error in if ((X - (j - 1) * dS) > f[i, j]) f[i, j] = (X - (j - 1) * dS) > : missing value where TRUE/FALSE needed > > On Windows both the stable 2.2.1 and 2.3.0rc gui versions will silently > produce incorrect data.Why do you file a bug report for a bug in your own code? The two loops can not be expected to give the same result since the term f[i,2] might differ when j is in 3:M. The platform difference could easily be due to floating point issues, e.g. the difference between divide by zero and divide by near-zero.> Nikolai. > > --------------020909040800030906040005 > Content-Type: text/plain; > name="broken.txt" > Content-Transfer-Encoding: 7bit > Content-Disposition: inline; > filename="broken.txt" > > findiff<-function(T,S_0,X,sigma,r,N,M) > { > f=array(,c(N+1,M+1)) > dT=T/N; > dS=2*S_0/M; > for (i in (1:(N+1))) > { > f[i,1]=X; > f[i,M+1]=0; > } > for (j in (1:(M+1))) > { > f[N+1,j]=max(X-(j-1)*dS,0); > } > a=0; > b=0; > c=0; > for (j in 1:M-1) > { > a[j+1]=.5*r*j*dT-.5*sigma^2*j^2*dT; > b[j+1]=1+sigma^2*j^2*dT+r*dT; > c[j+1]=-.5*r*j*dT-.5*sigma^2*j^2*dT; > } > i=N; > while (i>0) > { > d=0; > e=0; > d[1]=f[i+1,1]; > e[1]=0; > d[2]=0; > e[2]=1; > for (j in 2:M) > { > d[j+1]=(f[i+1,j]-a[j]*d[j-1]-b[j]*d[j])/c[j]; > e[j+1]=(-a[j]*e[j-1]-b[j]*e[j])/c[j]; > } > f[i,2]=(f[i,M+1]-d[M+1])/e[M+1]; > for (j in 2:M) > { > f[i,j]=d[j]+e[j]*f[i,2]; > if ((X-(j-1)*dS)>f[i,j]) f[i,j] = (X-(j-1)*dS); > } > i=i-1; > } > f; > } > > findiff(1,10,10,.3,.1,2,4) > > --------------020909040800030906040005 > Content-Type: text/plain; > name="working.txt" > Content-Transfer-Encoding: 7bit > Content-Disposition: inline; > filename="working.txt" > > findiff<-function(T,S_0,X,sigma,r,N,M) > { > f=array(,c(N+1,M+1)) > dT=T/N; > dS=2*S_0/M; > for (i in (1:(N+1))) > { > f[i,1]=X; > f[i,M+1]=0; > } > for (j in (1:(M+1))) > { > f[N+1,j]=max(X-(j-1)*dS,0); > } > a=0; > b=0; > c=0; > for (j in 1:M-1) > { > a[j+1]=.5*r*j*dT-.5*sigma^2*j^2*dT; > b[j+1]=1+sigma^2*j^2*dT+r*dT; > c[j+1]=-.5*r*j*dT-.5*sigma^2*j^2*dT; > } > i=N; > while (i>0) > { > d=0; > e=0; > d[1]=f[i+1,1]; > e[1]=0; > d[2]=0; > e[2]=1; > for (j in 2:M) > { > d[j+1]=(f[i+1,j]-a[j]*d[j-1]-b[j]*d[j])/c[j]; > e[j+1]=(-a[j]*e[j-1]-b[j]*e[j])/c[j]; > } > f[i,2]=(f[i,M+1]-d[M+1])/e[M+1]; > for (j in 2:M) > { > f[i,j]=d[j]+e[j]*f[i,2]; > } > for (j in 2:M) > { > if ((X-(j-1)*dS)>f[i,j]) f[i,j]=X-(j-1)*dS; > } > i=i-1; > } > f; > } > > findiff(1,10,10,.3,.1,2,4) > > --------------020909040800030906040005-- > > ______________________________________________ > R-devel at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel >-- O__ ---- Peter Dalgaard ?ster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
Hi Nicolai, 2006/4/20, N.Kalosha at math.ru.nl <N.Kalosha at math.ru.nl>:> This is a multi-part message in MIME format. > --------------020909040800030906040005 > Content-Type: text/plain; charset=KOI8-R; format=flowed > Content-Transfer-Encoding: 7bit > > Hi, > > I've attached two files with the sources for a function to implement the > finite difference method for solving a particular differential equation. > > One of them works correctly, the other gives wrong results or returns an > error, depending on the version of R. > > The difference between them is that in the 'broken' version in line 42 I > check if the items in the two-dimensional array are bigger than a > certain value, and in the working one I do it in a separate loop.diff working.r broken.r also reveals that your expression is differently bracketed broken.txt says in line 42: (...) = (X-(j-1)*dS); working.txt: (...) = X-(j-1)*dS; Don't know what the rationale of this difference really is since I'm no expert for finite difference methods. But your examples apparently don't only differ in terms of checking the size of your array-values. Does this matter? Thomas