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

[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.

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.

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