




Fortran Programming Language









unknown problem
Hi everybody I encounter with strange problem !!!! The below code scans plane and find the nearest value to zero . this works nice before that , but I there is no logical error with it !!!! Initial value of scan xaxis ( named xi , in program ) affected the answer !! If xi=0.5 the answer is different from when xi=1.5 . in both cases program find equal answers , but in last part when it searches for nearest value it makes a mistake . I don't it really doesn't depend on xi , but I can't find the source of error program test implicit none !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real,dimension(:,:),allocatable::A1(:,:) real, dimension(:,:),allocatable::A3(:,:) real,dimension(:),allocatable::A4(:) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real det ! determinant of jacobian real A ! ampliction real mag ! magnitude real u1,u2,u ! position of source real tE ! velocity of source . real b1 ! impcat parameter real theta ! the angle between sorce trace and horzintal axis real t0 !time of nearest ray to lens which centered on origin.in unit of day real pi real f, g real x,xi,xf,stepx,nx,ix real y,yi,yf,stepy,ny,iy real ti,tf,it,t,nt,stept real fmin,xans,yans,least integer i,j,n,m,k,p,ans,tot ! counter real m1,m2,m3,B,C,d,k1,k2 real x2,y2,x3,y3 !!!!!!!!!!!!! pi=4.0*atan(1.0) theta=pi/3 t0=0.2 tE=10 b1=0.5 m1=0.3 m2=0.3 m3=1(m1+m2) B=0 C=0 d =0.5 k1=0.5 k2=0.5 xi=0.5 xf=1.5 stepx=0.0001D0 yi=1. yf=1.5 stepy=0.0001D0 least=0.001D0 ti=0.75 tf=2. stept=0.02 i=0 ans=0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! allocate (A1(3,90000)) allocate (A3(3,90000)) allocate (A4(90)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!! t=0.35 u1=b1*sin(theta)+((tt0)/tE)*cos(theta) nx=(xfxi)/stepx ny=(yfyi)/stepy do ix=1,nx x=xi+(ix1)*stepx do iy=1,ny y=yi+(iy1)*stepy x2=xm1*B*x/(x**2+y**2) y2=ym1*B*y/(x**2+y**2) x3=x2m2*C*(x2d)/((x2d)**2+y2**2) y3=y2m2*C*y2/((x2d)**2+y2**2) f=abs(u1(xm1*x/(x**2+y**2)m2*(x2d)/((x2d)**2+y2**2)m3*(x3 k1)/((x3k1)**2+(y3k2)**2))) if (f.le.least) then u2=b1*cos(theta)((tt0)/tE)*sin(theta) g=abs(u2(ym1*y/(x**2+y**2)m2*y2/((x2d)**2+y2**2) m3*(y3k2)/((x3k1)**2+(y3k2)**2))) if (g.le.least)then i=i+1 A1(1,i)=x A1(2,i)=y A1(3,i)=f ! write (2,*),x,y,f end if ! if of g end if ! if of f end do ! do of x end do ! do of y ans=i n=0 do j=1,ans ! is the number of answers if (abs(A1(1,j)A1(1,j+1))>0.1.or.(abs(A1(2,j)A1(2,j+1))>0.1)) then n=n+1 A4(n)=j end if end do m=1 fmin=0.01D0 k=0 i=5 do j=m,A4(i) if (A1(3,j).lt.fmin)then write (*,*),j,fmin,A1(3,j) xans=A1(1,j) yans=A1(2,j) fmin=A1(3,j) end if end do ! do of j k=k+1 A3(1,k)=xans A3(2,k)=yans A3(3,k)=fmin write (*,*),"" write(*,*),xans,yans m=1+A4(i) end program best ,nakisa
The problem of detecting if lines intersect or not requires measuring a separation as being within a defined small amount. You have to subtract the variables, take the absolute value and test if this is within the acceptable distance. It is a common error to test if one variable is the same as another or to subtract and round truncate or convert to integers and test. There is always an exception that breaks the validity of the test.
On May 25, 11:51 am, nakisa <nakisa.noor@gmail.com> wrote: > Hi everybody > I encounter with strange problem !!!! > The below code scans plane and find the nearest value to zero . this > works nice before that , but I there is no logical error with > it !!!!
Oh yes there is. What exactly are you trying to minimize over a region in the XY plane? Is it f? > Initial value of scan xaxis ( named xi , in program ) affected the > answer !! > If xi=0.5 the answer is different from when xi=1.5 . in both > cases program find equal answers , but in last part when it searches > for nearest value it makes a mistake . > I don't it really doesn't depend on xi , but I can't find the source > of error
Well, if you are searching different regions of R2, you might get different answers. More to the point, I don't understand what you consider an "answer". First, a few stylistic comments on the program before looking at its "logic". > program test > implicit none > !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! > real,dimension(:,:),allocatable::A1(:,:) > real, dimension(:,:),allocatable::A3(:,:) > real,dimension(:),allocatable::A4(:) > !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
There is no apparent reason to make these arrays allocatable instead of ordinary arrays. > real det ! determinant of jacobian > real A ! ampliction > real mag ! magnitude
No need for stuff like this, but it does no real harm except lengthen the source listing.
> real u1,u2,u ! position of source > real tE ! velocity of source . > real b1 ! impcat parameter > real theta ! the angle between sorce trace and horzintal axis > real t0 !time of nearest ray to lens which centered on origin.in unit > of day > real pi > real f, g > real x,xi,xf,stepx,nx,ix > real y,yi,yf,stepy,ny,iy > real ti,tf,it,t,nt,stept > real fmin,xans,yans,least > integer i,j,n,m,k,p,ans,tot ! counter > real m1,m2,m3,B,C,d,k1,k2 > real x2,y2,x3,y3 > !!!!!!!!!!!!! > pi=4.0*atan(1.0) > theta=pi/3 > t0=0.2 > tE=10 > b1=0.5 > m1=0.3 > m2=0.3 > m3=1(m1+m2) > B=0 > C=0 > d =0.5 > k1=0.5 > k2=0.5 > xi=0.5 > xf=1.5 > stepx=0.0001D0
Assigning a double precision constant to a single precision real variable does not make the value which the variable contains more precise.
> yi=1. > yf=1.5 > stepy=0.0001D0 > least=0.001D0 > ti=0.75 > tf=2. > stept=0.02 > i=0 > ans=0 > !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! > allocate (A1(3,90000)) > allocate (A3(3,90000)) > allocate (A4(90)) > !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
See comment above. > t=0.35 > u1=b1*sin(theta)+((tt0)/tE)*cos(theta)
U2 is a loop invariant too. u2 =.......... belongs here instead of inside the loop. > nx=(xfxi)/stepx > ny=(yfyi)/stepy
nx, ny, ix and iy are REAL variables in this code. Using REAL instead of INTEGER variables in a DO loop is a relic of Fortran 77 code. Generally it is a "BAD IDEA". There is no reason not to use INTEGER, which is standard in Fortran 95 instead. > do ix=1,nx > x=xi+(ix1)*stepx > do iy=1,ny > y=yi+(iy1)*stepy > x2=xm1*B*x/(x**2+y**2) > y2=ym1*B*y/(x**2+y**2) > x3=x2m2*C*(x2d)/((x2d)**2+y2**2) > y3=y2m2*C*y2/((x2d)**2+y2**2)
Note that B and C are both zero, so these expressions are overly complex. Also note that x, x2, and x3 have the same value as do y, y2 and y3. > f=abs(u1(xm1*x/(x**2+y**2)m2*(x2d)/((x2d)**2+y2**2)m3*(x3 > k1)/((x3k1)**2+(y3k2)**2))) > if (f.le.least) then > u2=b1*cos(theta)((tt0)/tE)*sin(theta)
U2 takes on the same value each time it is calculated. This is called a "loop invariant". While leaving it here does not change the execution time of the program much, it really belongs outside the loop. > g=abs(u2(ym1*y/(x**2+y**2)m2*y2/((x2d)**2+y2**2) > m3*(y3k2)/((x3k1)**2+(y3k2)**2))) > if (g.le.least)then > i=i+1 > A1(1,i)=x > A1(2,i)=y > A1(3,i)=f > ! write (2,*),x,y,f
Uncommenting this statement and looking at the data file it generates would be instructive. > end if ! if of g > end if ! if of f > end do ! do of x
really do of y > end do ! do of y
really do of x > ans=i > n=0 > do j=1,ans ! is the number of answers > if (abs(A1(1,j)A1(1,j+1))>0.1.or.(abs(A1(2,j)A1(2,j+1))>0.1)) then > n=n+1 > A4(n)=j > end if > end do
Here you seem to be looking for the borders of regions where functions f and g are "small".
> m=1 > fmin=0.01D0 > k=0 > i=5 > do j=m,A4(i) > if (A1(3,j).lt.fmin)then > write (*,*),j,fmin,A1(3,j) > xans=A1(1,j) > yans=A1(2,j) > fmin=A1(3,j) > end if > end do ! do of j > k=k+1 > A3(1,k)=xans > A3(2,k)=yans > A3(3,k)=fmin > write (*,*),"" > write(*,*),xans,yans > m=1+A4(i) > end program
The rest of this code is obscure. Perhaps you need to move the end do lower in the listing, but I doubt that this will do what you want. Are you tying to find local minima in each "region" where abs(f) < tolerance and abs(g) < tolerance are both satisfied? HTH.  Elliot
hi . Dear Terence i need minimum values of abs(f) and abs(g) simulatiously . if i try abs(ffmin)<0.01 it means 0.01<ffmin<0.01 but it isnot thing i want. i search for the nearest value of zero of f and g . Dear Elliot at first ,thanks for your comments .i am new in fortran. i must find minimum value of " f" and "g "( both of them ).in other word, i must solve a system of 2 equations. surely , if i search different regions of R2 ,i may find different answers .but the problem is in for both case , where xi=0.5 and xi=1.5 , answer of first loop which scan R2 for minimum of f and g are the same . it is compeletly reasonable. but on selecting the least answer error arise ! i meas this program has error on third loop !! about B=0 and C=0, you are right . but this is first steps for trying this code .it must work for variouse values for B and C. you are right , i am trying to find local minima of f and g in each regions .the problem only take place on the last local , (i=5) . instead of "i=5" in original code was a loop .i remove some end part of this code , it is very long . the point which amaziang me is here . what is the relation between xi and xans when xi is 1 and xans is 1 !!!! i don't see it ! best ,nakisa
On May 26, 10:21 am, nakisa <nakisa.noor@gmail.com> wrote:
> hi . > Dear Terence > i need minimum values of abs(f) and abs(g) simulatiously . if i try > abs(ffmin)<0.01 it means > 0.01<ffmin<0.01 > but it isnot thing i want. i search for the nearest value of zero of > f and g . > Dear Elliot > at first ,thanks for your comments .i am new in fortran. > i must find minimum value of " f" and "g "( both of them ).in other > word, i must solve > a system of 2 equations. > surely , if i search different regions of R2 ,i may find different > answers .but the problem > is in for both case , where xi=0.5 and xi=1.5 , answer of first loop > which scan R2 for minimum > of f and g are the same . it is compeletly reasonable. but on > selecting the least answer error > arise ! i meas this program has error on third loop !! > about B=0 and C=0, you are right . but this is first steps for trying > this code .it must work > for variouse values for B and C. > you are right , i am trying to find local minima of f and g in each > regions .the problem only take place > on the last local , (i=5) . instead of "i=5" in original code was a > loop .i remove some end > part of this code , it is very long . > the point which amaziang me is here . what is the relation between xi > and xans when xi is 1 and > xans is 1 !!!! i don't see it ! > best ,nakisa
I made a number of changes to your original program. This is my corrected version. The changes are 1. change loop variables from REAL to INTEGER 2. save minimum values of BOTH f AND g  this required changing the size of A1 & A3 and adding the variable GMIN 3. calculation of u2 is moved out of the loop 4. finding "regions" accesses uninitialized data at (ans+1), so I set it to a "sentinal" value. 5. code restructured at end with nested loops For each region (I) search from M to A4(I). Before each search, set fmin and gmin to f and g of the first entry in that range. After each search, set M to A4(I)+1. program test implicit none !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real,dimension(:,:),allocatable::A1(:,:) real, dimension(:,:),allocatable::A3(:,:) real,dimension(:),allocatable::A4(:) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real det ! determinant of jacobian real A ! ampliction real mag ! magnitude real u1,u2,u ! position of source real tE ! velocity of source . real b1 ! impcat parameter real theta ! the angle between sorce trace and horzintal axis real t0 !time of nearest ray to lens which centered on origin.in unit of day real pi real f, g real x,xi,xf,stepx integer nx,ix ! was REAL real y,yi,yf,stepy integer ny,iy ! was REAL real ti,tf,it,t,nt,stept real fmin,xans,yans,least real gmin ! ADDED integer i,j,n,m,k,p,ans,tot ! counter real m1,m2,m3,B,C,d,k1,k2 real x2,y2,x3,y3 !!!!!!!!!!!!! pi=4.0*atan(1.0) theta=pi/3 t0=0.2 tE=10 b1=0.5 m1=0.3 m2=0.3 m3=1(m1+m2) B=0 C=0 d =0.5 k1=0.5 k2=0.5 xi=0.5 xf=1.5 stepx=0.0001D0 yi=1. yf=1.5 stepy=0.0001D0 least=0.001D0 ti=0.75 tf=2. stept=0.02 i=0 ans=0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! allocate (A1(4,90000)) ! was 3,90000 allocate (A3(4,90000)) ! was 3,90000 allocate (A4(90)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!! t=0.35 u1=b1*sin(theta)+((tt0)/tE)*cos(theta) u2=b1*cos(theta)((tt0)/tE)*sin(theta) ! moved outside of loop nx=(xfxi)/stepx ny=(yfyi)/stepy do ix=1,nx x=xi+(ix1)*stepx do iy=1,ny y=yi+(iy1)*stepy x2=xm1*B*x/(x**2+y**2) y2=ym1*B*y/(x**2+y**2) x3=x2m2*C*(x2d)/((x2d)**2+y2**2) y3=y2m2*C*y2/((x2d)**2+y2**2) f=abs(u1(xm1*x/(x**2+y**2)m2*(x2d)/((x2d)**2+y2**2)m3*(x3 k1)/ ((x3k1)**2+(y3k2)**2))) if (f.le.least) then ! u2=b1*cos(theta)((tt0)/tE)*sin(theta) g=abs(u2(ym1*y/(x**2+y**2)m2*y2/((x2d)**2+y2**2)m3*(y3k2)/ ((x3 k1)**2+(y3k2)**2))) if (g.le.least)then i=i+1 A1(1,i)=x A1(2,i)=y A1(3,i)=f A1(4,i)=g ! save g as well ! write (2,*) x,y,f,g end if ! if of g end if ! if of f end do ! do of x end do ! do of y A1(1,i+1)=999 ! end of data + 1 A1(2,i+1)=999 ans=i n=0 do j=1,ans ! is the number of answers if (abs(A1(1,j)A1(1,j+1))>0.1.or.(abs(A1(2,j)A1(2,j+1))>0.1)) then n=n+1 A4(n)=j end if end do ! restructured from here on m=1 do i=1,n fmin=A1(3,m) ! set to first gmin=A1(4,m) ! set to first do j=m,A4(i) if (A1(3,j).lt.fmin) then if (A1(4,j).lt.gmin) then ! write (*,*) i,j,A1(1,j),A1(2,j),A1(3,j),A1(4,j),fmin,gmin xans=A1(1,j) yans=A1(2,j) fmin=A1(3,j) gmin=A1(4,j) ! save g as well end if ! of gmin end if ! of fmin end do ! do of j A3(1,i)=xans A3(2,i)=yans A3(3,i)=fmin A3(4,i)=gmin ! write (*,*),"" write (*,*) xans,yans m=1+A4(i) end do ! do of i end program With these changes, the output was 0.3578000 0.4837000 0.1571000 1.220300 0.2272000 4.0299974E02 0.4577000 0.2438000 1.105400 0.9751999 As expected, there are 5 regions, with one local minimum for each of them. This output was the same even if xi was set to 1.5 instead of .5. HTH  Elliot
e p chandler wrote:
... > real,dimension(:,:),allocatable::A1(:,:) > real, dimension(:,:),allocatable::A3(:,:) > real,dimension(:),allocatable::A4(:)
Well, it doesn't mean anything different, but why the redundant shape specifications of the variables? The above means the same as either: real,dimension(:,:),allocatable::A1 real, dimension(:,:),allocatable::A3 real,dimension(:),allocatable::A4 or: real,allocatable::A1(:,:) real,allocatable::A3(:,:) real,allocatable::A4(:) I prefer the latter. In fact, I think the DIMENSION keyword should be removed from the language. The supposed advantage is that it allows you to simultaneously declare several arrays the same shape (which you don't even do here). If you are going to use a separate declaration statement for each variable, it serves no purpose at all.  J. Giles "I conclude that there are two ways of constructing a software design: One way is to make it so simple that there are obviously no deficiencies and the other way is to make it so complicated that there are no obvious deficiencies."  C. A. R. Hoare
On May 26, 3:32 pm, "James Giles" <jamesgi@worldnet.att.net> wrote:
> e p chandler wrote: > ... > > real,dimension(:,:),allocatable::A1(:,:) > > real, dimension(:,:),allocatable::A3(:,:) > > real,dimension(:),allocatable::A4(:) > Well, it doesn't mean anything different, but why the redundant > shape specifications of the variables? The above means the > same as either: > real,dimension(:,:),allocatable::A1 > real, dimension(:,:),allocatable::A3 > real,dimension(:),allocatable::A4 > or: > real,allocatable::A1(:,:) > real,allocatable::A3(:,:) > real,allocatable::A4(:) > I prefer the latter. In fact, I think the DIMENSION keyword > should be removed from the language. The supposed advantage > is that it allows you to simultaneously declare several arrays > the same shape (which you don't even do here). If you are > going to use a separate declaration statement for each variable, > it serves no purpose at all. >  > J. Giles
Ouch, it's not my code! I was reworking the OP's program. Ignorance, yes. I don't know enough Fortran 95 to use ALLOCATABLE in my own programs, so I did not see the redundancy. I would not have used ALLOCATABLE in this program, but after making a number of changes to get it to work and compile, I gave up "tweaking".  Elliot
"e p chandler" <e@juno.com> wrote in message news:1180214913.258361.128580@o5g2000hsb.googlegroups.com... > On May 26, 3:32 pm, "James Giles" <jamesgi @worldnet.att.net> wrote: ... >> I prefer the latter. In fact, I think the DIMENSION keyword >> should be removed from the language. The supposed advantage >> is that it allows you to simultaneously declare several arrays >> the same shape (which you don't even do here). If you are >> going to use a separate declaration statement for each variable, >> it serves no purpose at all. ... > Ouch, it's not my code! I was reworking the OP's program. Ignorance, > yes. I don't know enough Fortran 95 to use ALLOCATABLE in my own > programs, so I did not see the redundancy. I would not have used > ALLOCATABLE in this program, but after making a number of changes to > get it to work and compile, I gave up "tweaking". I wasn't criticizing you particularly, just the code. I would have responded to the OP if I'd seen that first. Sorry to make it seem directed at you.  J. Giles "I conclude that there are two ways of constructing a software design: One way is to make it so simple that there are obviously no deficiencies and the other way is to make it so complicated that there are no obvious deficiencies."  C. A. R. Hoare
On May 26, 10:21 am, nakisa <nakisa.noor@gmail.com> wrote:
> hi . > Dear Terence > i need minimum values of abs(f) and abs(g) simulatiously . if i try > abs(ffmin)<0.01 it means > 0.01<ffmin<0.01 > but it isnot thing i want. i search for the nearest value of zero of > f and g . > Dear Elliot > at first ,thanks for your comments .i am new in fortran. > i must find minimum value of " f" and "g "( both of them ).in other > word, i must solve > a system of 2 equations. > surely , if i search different regions of R2 ,i may find different > answers .but the problem > is in for both case , where xi=0.5 and xi=1.5 , answer of first loop > which scan R2 for minimum > of f and g are the same . it is compeletly reasonable. but on > selecting the least answer error > arise ! i meas this program has error on third loop !! > about B=0 and C=0, you are right . but this is first steps for trying > this code .it must work > for variouse values for B and C. > you are right , i am trying to find local minima of f and g in each > regions .the problem only take place > on the last local , (i=5) . instead of "i=5" in original code was a > loop .i remove some end > part of this code , it is very long . > the point which amaziang me is here . what is the relation between xi > and xans when xi is 1 and > xans is 1 !!!! i don't see it ! > best ,nakisa
OK. Refer to other posts in this thread, including those by James Giles for changes that you should make to your program. I'm going to try and answer your original question now that I understand the problem better. Your X values are going to differ if you change the starting point from .5 to 1.5. This is a consequence of the properties of floating point arithmetic, which almost all of the time does not deal with exact numbers. Your step size of .0001 does not have an exact binary representation inside the computer. Instead it is stored as an approximate value. So when you compute x = xi + (ix1)*stepx each x has a small error associated with it because stepx is really approximate. Note that when you increase your search range by lowering the start of the search from .5 to 1.5, you increase the number of steps taken. IX is LARGER near the end of the range then it was before. So the ERROR in x increases. The result is that your x values are actually different in the two runs. HTH Back to the problem your program is trying to solve. It might be better to use this type of program with a smaller step size to find regions where F and G are a minimum. Then write a program which uses a two dimensional version of Newton's method to find ROOTS in (x,y) given the solution the particular regions you are interested in as an initial quess.  elliot  elliot
dear elliot thanks alot of your helps best,nakisa
hi i have another question . when the program started, it scans R2 with 0.01 precision on x and y, but the answers are real numbers with 4 digit after point. these answers will be used in other parts of program and errors by new digits are unavoidable !! i think may be there is a way to cut unwanted digits , such as in various formats in printing data . is there any way ? best ,nakisa
Re: solving 2 nonlinear equations in two unknowns On May 27, 11:06 pm, nakisa <nakisa.noor@gmail.com> wrote: > dear elliot > thanks alot of your helps > best,nakisa
You are welcome. recap: The original program searched over a range in (x,y) to find points where functions f and g were minimal. From these, it isolated regions. In each region, it found (x,y) with the smallest values of both f and g. An alternate approach is to apply a two dimensional version of Newton's method. I simplified the functions f and g by removing abs(), and setting B & C to 0. We really are looking for approximate roots (x,y) of f & g, so removing the abs() is justified. This leaves: f=u1(xm1*x/(x**2+y**2)m2*(xd)/((xd)**2+y**2) & m3*(xk1)/((xk1)**2+(yk2)**2)) g=u2(ym1*y/(x**2+y**2)m2*y/((xd)**2+y**2) & m3*(yk2)/((xk1)**2+(yk2)**2)) This is an iterative process involving an initial guess, functions f & g, and the inverse of the Jacobean matrix. (partial derivatives of f & g with respect to x & y). At each step, I calculated the Jacobean and took its inverse using a well known formula for the inverse of a 2x2 matrix. [Flip the matrix about the minor diagonal, negate the entries on the minor diagonal and divide by the determinant of the original matrix.] Repeat until the "difference" between old (x,y) and new (x,y) is "small" enough. Two methods both gave similar results: a. approximate partial derivatives  (f(x+h,y)f(x,y)) / h, for suitable small value of h, etc. b. "exact" partial derivatives  from elementary calculus [messy but not difficult] Here is a table of results 0.3578000 0.4837000 (.3,.4) 0.3581022 0.4834902 3.7145835E08 3.2571472E07 0.3581020 0.4834905 1.3417747E07 2.8465786E07 0.1571000 1.220300 (.2,1) 0.1562073 1.219877 2.3027271E08 9.9773780E09 0.1562074 1.219878 2.8230591E08 1.9128676E07 0.2272000 4.0299974E02 (.2,0) 0.2272325 4.0326063E02 1.8983251E08 2.9303777E07 0.2272325 4.0326044E02 1.2422674E07 1.2010982E07 0.4577000 0.2438000 (.45,.25) 0.4576834 0.2438498 2.8040859E07 6.1542141E07 0.4576833 0.2438498 8.1046187E07 6.9876904E07 1.105400 0.9751999 (1.1,.975) 0.1562070 1.219877 2.4550855E09 1.9470619E08 0.1562067 1.219877 7.8551810E10 8.8881983E08 Each group corresponds to a "region" found by my original program. The first line lists its (x,y) and my initial guess. Lines 2 and 3 are the results from approximate partial derivatives and "exact" partial derivatives  listing x, y, f and g. Note that I could not get the algorithm to converge to the (x,y) in region 5. As you see the program drifted into a different region. (True for various initial (x,y).) Here is part of the 2nd program. It uses statement functions. It's ugly. [snip] fnf(x,y)=u1(xm1*x/(x**2+y**2)m2*(xd)/((xd)**2+y**2) & m3*(xk1)/((xk1)**2+(yk2)**2)) fng(x,y)=u2(ym1*y/(x**2+y**2)m2*y/((xd)**2+y**2) & m3*(yk2)/((xk1)**2+(yk2)**2)) [snip] i=0 do f=fnf(x,y) g=fng(x,y) fx=fnf(x+h,y) fy=fnf(x,y+h) gx=fng(x+h,y) gy=fng(x,y+h) q=(fxf)*(gyg)  (fyf)*(gxg) if ( abs(q) < qtol ) then print *,'determinant ~ zero' exit end if x1=x  h/q * ( f*(gyg)  g*(fyf) ) ! X* = X  Jinv F y1=y  h/q * (f*(gxg) + g*(fxf) ) i=i+1 if ( sqrt( (x1x)**2 + (y1y)**2 ) < tol ) then [snip] Here is the whole third program: implicit none real u1,u2,tE,b1,theta,t0,pi,t real m1,m2,m3,d,k1,k2 real f,g,fx,gx,fy,gy real x,x1 real y,y1 real tol,qtol real q integer i,max pi=4.0*atan(1.0) theta=pi/3 t0=0.2 tE=10 b1=0.5 t=0.35 u1=b1*sin(theta)+((tt0)/tE)*cos(theta) u2=b1*cos(theta)((tt0)/tE)*sin(theta) m1=0.3 m2=0.3 m3=1(m1+m2) d=0.5 k1=0.5 k2=0.5 tol=1e6 max=30 qtol=1e15 print *,'x0,y0:' read *,x,y i=0 do f=u1(xm1*x/(x**2+y**2)m2*(xd)/((xd)**2+y**2) & m3*(xk1)/((xk1)**2+(yk2)**2)) g=u2(ym1*y/(x**2+y**2)m2*y/((xd)**2+y**2) & m3*(yk2)/((xk1)**2+(yk2)**2)) fx= 1 & m1*(x**2y**2)/(x**2+y**2)**2 & m2*((xd)**2y**2)/((xd)**2+y**2)**2 & m3*((xk1)**2(yk2)**2)/((xk1)**2+(yk2)**2)**2 fy= 2*m1*x*y/(x**2+y**2)**2 & 2*m2*(xd)*y/((xd)**2+y**2)**2 & 2*m3*(xk1)*(yk2)/((xk1)**2+(yk2)**2)**2 gx= 2*m1*x*y/(x**2+y**2)**2 & 2*m2*(xd)*y/((xd)**2+y**2)**2 & 2*m3*(xk1)*(yk2)/((xk1)**2+(yk2)**2)**2 gy= 1 & +m1*(x**2y**2)/(x**2+y**2)**2 & +m2*((xd)**2y**2)/((xd)**2+y**2)**2 & +m3*((xk1)**2(yk2)**2)/((xk1)**2+(yk2)**2)**2 q=fx*gy  fy*gx ! determinant of Jacobean if ( abs(q) < qtol ) then print *,'determinant ~ zero' exit end if x1=x  1/q * ( f*gy  g*fy ) ! X* = X  Jinv F y1=y  1/q * (f*gx + g*fx ) i=i+1 print *,i,x,y,x1,y1 if ( sqrt( (x1x)**2 + (y1y)**2 ) < tol ) then print *,'x,y=',x,y print *,'f,g=',f,g exit end if x=x1 y=y1 if (i == max) then print *,'convergence failed at step ',max exit end if end do end Note that the equations in your original program are more complex. :).  Elliot
nakisa <nakisa.noor @gmail.com> wrote: > but the answers > are real numbers with 4 digit after point. No, they are not. You really, really need to learn some of the basics of computer floating point (not specific to Fortran). Real numbers are not decimal at all. It is pretty much meaningless to talk about how many decimal digits they have. They don't have decimal digits because they aren't in decimal. They are invariably in binary. Anything that you see in a decimal form is the result of converting the binary value to an approximate decimal equivalent for display purposes. It is incredibly, incredibly important to understand this. It is at the core of all current floating point arithmetic. (Some people here might point out the theoretical possibility of decimal floating point, and even that some machine might hav ethe theoretical capability for it. But that is pretty much irrelevant for your purposes. Anything you will see will be binary floating point.)  Richard Maine  Good judgement comes from experience; email: last name at domain . net  experience comes from bad judgement. domain: summertriangle   Mark Twain
hi Elliot thanks alot nakisa





