C*************************************************************** C C MAXFUN C ============ C C Subrutina de calcul a maximului unei functii neliniare C de o variabila. C Variabila nu este supusa nici unei restrictii. C Nu se utilizeaza derivata functiei pentru care se cauta C punctul de maxim. C C max(f(x)) C C Subrutina implementeaza metoda de cautare directa Powell. C C Secventa de apel: C call maxfun(eps,pas,x1) C unde: C eps = acuratetea de calcul a solutiei. C pas = lungimea pasului de cautare. C x1 = punctul initial (de unde incepe cautarea maximului). C C Subrutine apelate: C Subroutine func(x,y) C x punctul in care se calculeaza valoarea functiei C y = f(x) valoarea functiei in punctul x. C C Neculai Andrei C 1980, ICI - Bucuresti C*************************************************************** subroutine maxfun(eps,pas,x1) ** double precision t1(3), t2(3), vf(3) double precision eps,pas double precision x1,y1, x2,y2, x3,y3, x4,y4 double precision a,b,c,d,e,f,p,r,xe,xmax,vfxe,vfmax logical iprint * Open the output file open(unit=2,file='maxfun.out',status='unknown') iprint=.true. call func(x1,y1) if(iprint) then write(2,100) 100 format(/,24x,1hX,20x,4hF(X),17x,5hPASUL) write(2,101) x1,y1,pas 101 format(15x,e20.13,2x,e20.13,2x,e20.13) end if * Incrementarea variabilei x2=x1+pas ii=0 1 call func(x2,y2) if(iprint) write(2,101) x2,y2,pas * Compararea valorilor functiei n cele doua puncte if(y2-y1) 3,6,6 3 if(ii-1) 4,4,5 4 ii=ii+1 * Schimbarea directiei de cautare pas=-pas x2=x1+pas go to 1 5 go to 11 6 continue x3=x2+pas go to 8 7 x3=x4 8 call func(x3,y3) if(iprint) write(2,101) x3,y3,pas * Dublarea pasului pas=pas*2.d0 x4=x3+pas call func(x4,y4) if(y4-y3) 9,10,10 9 go to 12 10 go to 7 * Evaluarea a trei puncte in jurul maximului 11 t1(1)=x1 t1(2)=x1+pas/2.d0 t1(3)=x2 go to 13 12 t1(1)=x3 t1(2)=x3+pas/2.d0 t1(3)=x4 go to 13 * Evaluarea functiei in aceste puncte 13 do 14 i=1,3 call func(t1(i),vf(i)) 14 continue * Interpolarea patratica utilizand aceste puncte if(iprint) write(2,102) 102 format(/,24x,2hX*,19x,5hF(X*)) 15 continue a=t1(2)-t1(3) b=t1(3)-t1(1) c=t1(1)-t1(2) d=t1(2)**2-t1(3)**2 e=t1(3)**2-t1(1)**2 f=t1(1)**2-t1(2)**2 p=0.5d0*(d*vf(1)+e*vf(2)+f*vf(3)) r= a*vf(1)+b*vf(2)+c*vf(3) xe=p/r call func(xe,vfxe) if(iprint) write(2,103)xe,vfxe 103 format(15x,e20.13,2x,e20.13) do 18 j=1,3 t2(j)=dabs(t1(j)-xe) if(t2(j)-eps) 16,16,17 16 xmax=t1(j) go to 25 17 continue 18 continue * Determinarea si inlocuirea punctului de valoare maxima if(vf(2)-vf(1)) 20,19,19 19 if(vf(3)-vf(1)) 20,21,21 20 if(vf(3)-vf(2)) 23,22,22 21 jk=1 go to 24 22 jk=2 go to 24 23 jk=3 24 t1(jk)=xe vf(jk)=vfxe go to 15 25 call func(xmax,vfmax) write(2,104) xmax 104 format(/15x,'Punctul de maxim determinat. xmax=',2x,e20.13) write(2,105) vfmax 105 format(/15x,'Valoarea functiei in punctul xmax=',2x,e20.13) write(2,106) eps 106 format(/15x,'Acuratetea de calcul a solutiei ',2x,e20.13) return end **************************************************************** subroutine func(x,y) double precision x,y * Exemplul 1 c1 y = -x*x + 6.d0*x - 2.d0 c1 y =-y * Exemplul 2 c2 y = x**5 - 2.d0*x**3 +10.d0*dsin(5.d0*x) c2 y =-y * Exemplul 3 c3 y = 1.d0 - 10.d0*x + 0.01d0*exp(x) c3 y =-y * Exemplul 4 c4 y=-1.d0+2.d0/(x**2+1.d0) * Exemplul 5 c y=-x/(x*x+1.d0) * Exemplul 6 c6 y=5.d0*x**5-4.d0*x**4+400.d0*x*dsin(4.d0*x-4.d0) c6 y=-y * Exemplul 7 c7 y=x**4-12.d0*x**3+15.d0*x**2+56.d0*x-60.d0 * Exemplul 8 y=-x**4+12.d0*x**3-47.d0*x**2+60.d0*x return end **************************************************************** * Main Program for MAXFUN subroutine double precision eps,pas,x1 eps= 0.00001d0 pas= 0.001d0 x1 = -1.5d0 call maxfun(eps,pas,x1) stop end ****************************************************** Last line