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

From: Peter Dalgaard <p.dalgaard_at_biostat.ku.dk>
Date: Thu 20 Apr 2006 - 10:50:43 GMT

N.Kalosha@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@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@biostat.ku.dk)                  FAX: (+45) 35327907

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

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