SR1校正的拟牛顿法

SR1校正的拟牛顿法 0

!!!本程序适用于求解形如f(x)=1/2*x'Ax+bx+c二次函数的稳定点； !!!输入函数信息，输出函数的稳定点及迭代次数； !!!iter整型变量，存放迭代次数；

!!!x为n维变量，初始值由用户输入；gradt实型变量，存放函数梯度; !!!dir实型变量，存放搜索方向； program main

real,dimension(:),allocatable::x,gradt,dir,b ,s,y,p ,gradt1,x1 real,dimension(:,:),allocatable::hessin ,H ,G real::x0,tol integer::n ,iter,i,j print*,'请输入变量的维数' read*,n

allocate(x(n),gradt(n),dir(n),b(n),s(n),y(n),p(n),gradt1(n),x1(n)) allocate(hessin(n,n),H(n,n),G(n,n)) print*,'请输入初始向量x' read*,x

print*,'请输入hessin矩阵' read*,hessin print*,'请输入矩阵b' read*,b iter=0 tol=0.000001

do i=1,n do j=1,n if (i==j)then H(i,j)=1 else H(i,j)=0 endif enddo enddo

100 gradt=matmul(hessin,x)+b

if(sqrt(dot_product(gradt,gradt))dir=-matmul(H,gradt)

x0=golden(x,dir,hessin,b) x1=x+x0*dir

gradt1=matmul(hessin,x1)+b s=x1-x y=gradt1-gradt p=s-matmul(H,y) call vectorm(p,G)

H=H+1/dot_product(p,y)*G x=x1 iter=iter+1 if(iter>10*n)then print*,\"out\" goto 101 endif

print*,\"第\次运行结果为\方向为\步长\print*,x,\"f(x)=\ goto 100 contains

!!!子程序，返回函数值 function f(x,A,b) result(f_result) real,dimension(:),intent(in)::x,b real,dimension(:,:),intent(in)::A real::f_result

f_result=0.5*dot_product(matmul(x,A),x)+dot_product(b,x) end function f

!!!子程序，矩阵与向量相乘 subroutine vectorm(p,G) real,dimension(:),intent(in)::p real,dimension(:,:),intent(out)::G n=size(p) do i=1,n !do j=1,n G(i,:)=p(i)*p !enddo enddo

end subroutine

!!!精确线搜索0.618法子程序 ，返回步长； function golden(x,d,A,b) result(golden_n) real::golden_n real::x0

real,dimension(:),intent(in)::x,d real,dimension(:),intent(in)::b

real,dimension(:,:),intent(in)::A real::x1,x2,a1,b1,f0,f1,f2,r,tol,dx parameter(r=0.618) tol=0.0001 dx=0.1 x0=1 x1=x0+dx f0=f(x+x0*d,A,b) f1=f(x+x1*d,A,b) if(f0=f1)then b1=x2 a1=x0 else x0=x1 x1=x2 f0=f1 f1=f2 goto 2 endif endif

x1=a1+(1-r)*(b1-a1) x2=a1+r*(b1-a1) f1=f(x+x1*d,A,b) f2=f(x+x2*d,A,b) 3 if(abs(b1-a1)<=tol)then x0=(a1+b1)/2

else

if(f1>f2)then a1=x1 x1=x2 f1=f2

x2=a1+r*(b1-a1) f2=f(x+x2*d,A,b) goto 3 else b1=x2 x2=x1 f2=f1

x1=a1+(1-r)*(b1-a1) f1=f(x+x1*d,A,b) goto 3 endif endif golden_n=x0 end function golden

101 end

本算法由Fortran 90语言编写，在Visual Fortran 5上编译通过，本程序由沙沙提供！