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
