|
|
 |
 |
 |
 |
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(:)
> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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
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))
> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
> 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
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)
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)
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
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
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
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(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:
> 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
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:
> 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
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".
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(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
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.
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