Home     |     .Net Programming    |     cSharp Home    |     Sql Server Home    |     Javascript / Client Side Development     |     Ajax Programming

Ruby on Rails Development     |     Perl Programming     |     C Programming Language     |     C++ Programming     |     IT Jobs

Python Programming Language     |     Laptop Suggestions?    |     TCL Scripting     |     Fortran Programming     |     Scheme Programming Language


 
 
Cervo Technologies
The Right Source to Outsource

MS Dynamics CRM 3.0

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 x-axis ( 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)+((t-t0)/tE)*cos(theta)

   nx=(xf-xi)/stepx
   ny=(yf-yi)/stepy

       do ix=1,nx
             x=xi+(ix-1)*stepx

          do iy=1,ny
            y=yi+(iy-1)*stepy

                x2=x-m1*B*x/(x**2+y**2)
            y2=y-m1*B*y/(x**2+y**2)

            x3=x2-m2*C*(x2-d)/((x2-d)**2+y2**2)
            y3=y2-m2*C*y2/((x2-d)**2+y2**2)

                    f=abs(u1-(x-m1*x/(x**2+y**2)-m2*(x2-d)/((x2-d)**2+y2**2)-m3*(x3-
k1)/((x3-k1)**2+(y3-k2)**2)))

              if (f.le.least) then

                 u2=b1*cos(theta)-((t-t0)/tE)*sin(theta)
                 g=abs(u2-(y-m1*y/(x**2+y**2)-m2*y2/((x2-d)**2+y2**2)-
m3*(y3-k2)/((x3-k1)**2+(y3-k2)**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 X-Y plane? Is it f?

> Initial value of scan x-axis ( 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.

Assigning a double precision constant to a single precision real
variable does not make the value which the variable contains more
precise.

See comment above.

> t=-0.35
> u1=b1*sin(theta)+((t-t0)/tE)*cos(theta)

U2 is a loop invariant too.
u2 =..........
belongs here instead of inside the loop.

>    nx=(xf-xi)/stepx
>    ny=(yf-yi)/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+(ix-1)*stepx

>           do iy=1,ny
>             y=yi+(iy-1)*stepy

>                 x2=x-m1*B*x/(x**2+y**2)
>             y2=y-m1*B*y/(x**2+y**2)

>             x3=x2-m2*C*(x2-d)/((x2-d)**2+y2**2)
>             y3=y2-m2*C*y2/((x2-d)**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-(x-m1*x/(x**2+y**2)-m2*(x2-d)/((x2-d)**2+y2**2)-m3*(x3-
> k1)/((x3-k1)**2+(y3-k2)**2)))

>               if (f.le.least) then

>                  u2=b1*cos(theta)-((t-t0)/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-(y-m1*y/(x**2+y**2)-m2*y2/((x2-d)**2+y2**2)-
> m3*(y3-k2)/((x3-k1)**2+(y3-k2)**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

Un-commenting 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".

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?

> best ,nakisa

HTH.

-- Elliot

hi .

Dear Terence
i need minimum values of abs(f) and abs(g) simulatiously . if i try
abs(f-fmin)<0.01 it means
-0.01<f-fmin<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:

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 un-initialized 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)+((t-t0)/tE)*cos(theta)
u2=b1*cos(theta)-((t-t0)/tE)*sin(theta) ! moved outside of loop

nx=(xf-xi)/stepx
ny=(yf-yi)/stepy

do ix=1,nx
  x=xi+(ix-1)*stepx

  do iy=1,ny
    y=yi+(iy-1)*stepy

    x2=x-m1*B*x/(x**2+y**2)
    y2=y-m1*B*y/(x**2+y**2)

    x3=x2-m2*C*(x2-d)/((x2-d)**2+y2**2)
    y3=y2-m2*C*y2/((x2-d)**2+y2**2)

    f=abs(u1-(x-m1*x/(x**2+y**2)-m2*(x2-d)/((x2-d)**2+y2**2)-m3*(x3-
k1)/

((x3-k1)**2+(y3-k2)**2)))
    if (f.le.least) then
!     u2=b1*cos(theta)-((t-t0)/tE)*sin(theta)
      g=abs(u2-(y-m1*y/(x**2+y**2)-m2*y2/((x2-d)**2+y2**2)-m3*(y3-k2)/
((x3-

k1)**2+(y3-k2)**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.0299974E-02
  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:

Ouch, it's not my code! I was re-working 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 re-working 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:

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 + (ix-1)*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 non-linear 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-(x-m1*x/(x**2+y**2)-m2*(x-d)/((x-d)**2+y**2) &
    -m3*(x-k1)/((x-k1)**2+(y-k2)**2))

  g=u2-(y-m1*y/(x**2+y**2)-m2*y/((x-d)**2+y**2) &
    -m3*(y-k2)/((x-k1)**2+(y-k2)**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.7145835E-08  3.2571472E-07
-0.3581020  -0.4834905       1.3417747E-07  2.8465786E-07

0.1571000   1.220300         (.2,1)
0.1562073   1.219877        -2.3027271E-08 -9.9773780E-09
0.1562074   1.219878         2.8230591E-08 -1.9128676E-07

0.2272000   4.0299974E-02    (.2,0)
0.2272325   4.0326063E-02   -1.8983251E-08  2.9303777E-07
0.2272325   4.0326044E-02   -1.2422674E-07  1.2010982E-07

0.4577000   0.2438000        (.45,.25)
0.4576834   0.2438498       -2.8040859E-07 -6.1542141E-07
0.4576833   0.2438498       -8.1046187E-07 -6.9876904E-07

1.105400    0.9751999        (1.1,.975)
0.1562070   1.219877        -2.4550855E-09  1.9470619E-08
0.1562067   1.219877        -7.8551810E-10  8.8881983E-08

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-(x-m1*x/(x**2+y**2)-m2*(x-d)/((x-d)**2+y**2) &
      -m3*(x-k1)/((x-k1)**2+(y-k2)**2))

fng(x,y)=u2-(y-m1*y/(x**2+y**2)-m2*y/((x-d)**2+y**2) &
      -m3*(y-k2)/((x-k1)**2+(y-k2)**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=(fx-f)*(gy-g) - (fy-f)*(gx-g)
  if ( abs(q) < qtol ) then
    print *,'determinant ~ zero'
    exit
  end if
  x1=x - h/q * ( f*(gy-g) - g*(fy-f) ) ! X* = X - Jinv F
  y1=y - h/q * (-f*(gx-g) + g*(fx-f) )
  i=i+1
  if ( sqrt( (x1-x)**2 + (y1-y)**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)+((t-t0)/tE)*cos(theta)
u2=b1*cos(theta)-((t-t0)/tE)*sin(theta)

m1=0.3
m2=0.3
m3=1-(m1+m2)
d=0.5
k1=0.5
k2=0.5

tol=1e-6
max=30
qtol=1e-15

print *,'x0,y0:'
read *,x,y

i=0
do
  f=u1-(x-m1*x/(x**2+y**2)-m2*(x-d)/((x-d)**2+y**2) &
    -m3*(x-k1)/((x-k1)**2+(y-k2)**2))

  g=u2-(y-m1*y/(x**2+y**2)-m2*y/((x-d)**2+y**2) &
    -m3*(y-k2)/((x-k1)**2+(y-k2)**2))

  fx= -1 &
    -m1*(x**2-y**2)/(x**2+y**2)**2 &
    -m2*((x-d)**2-y**2)/((x-d)**2+y**2)**2 &
    -m3*((x-k1)**2-(y-k2)**2)/((x-k1)**2+(y-k2)**2)**2

  fy= -2*m1*x*y/(x**2+y**2)**2 &
    -2*m2*(x-d)*y/((x-d)**2+y**2)**2 &
    -2*m3*(x-k1)*(y-k2)/((x-k1)**2+(y-k2)**2)**2

  gx= -2*m1*x*y/(x**2+y**2)**2 &
         -2*m2*(x-d)*y/((x-d)**2+y**2)**2 &
         -2*m3*(x-k1)*(y-k2)/((x-k1)**2+(y-k2)**2)**2

  gy= -1 &
    +m1*(x**2-y**2)/(x**2+y**2)**2 &
    +m2*((x-d)**2-y**2)/((x-d)**2+y**2)**2 &
    +m3*((x-k1)**2-(y-k2)**2)/((x-k1)**2+(y-k2)**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( (x1-x)**2 + (y1-y)**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
Add to del.icio.us | Digg this | Stumble it | Powered by Megasolutions Inc