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

From: <N.Kalosha_at_math.ru.nl>
Date: Thu 20 Apr 2006 - 10:26:10 GMT


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



R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Thu Apr 20 20:33:30 2006

This archive was generated by hypermail 2.1.8 : Thu 20 Apr 2006 - 12:17:46 GMT