[Rd] bug: code not working as expected (PR#8783)

[Peter Dalgaard]

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