[EMAIL PROTECTED] 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
~~~~~~~~~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
