拟牛顿BFS方法求解复杂四元非线性超函数方程组问题

无奈的小强

[Fortran] 纯文本查看 复制代码

module inv_mat

contains
subroutine inv(A,invA,N)
implicit real*8(a-z)
integer::n
integer::i
real*8::A(n,n),invA(n,n),E(n,n)
E=0
!设置E为单位矩阵
do i=1,n
   E(i,i)=1
end do
call mateq(A,E,invA,N,N)
end subroutine inv
subroutine mateq(A,B,X,N,M)
implicit real*8(a-z)
integer::N,M,i
real*8::A(N,N),B(N,M),X(N,M)
real*8::btemp(N),xtemp(N)
do i=1,M    
    btemp=B(:,i)
    call elgauss(A,btemp,xtemp,N)   
    X(:,i)=xtemp
end do
end subroutine mateq
subroutine elgauss(A,b,x,N)
implicit real*8(a-z)
integer::i,k,N
integer::id_max  !主元素标号
real*8::A(N,N),b(N),x(N)
real*8::Aup(N,N),bup(N)
!Ab为增广矩阵  [Ab]
real*8::Ab(N,N+1)
real*8::vtemp1(N+1),vtemp2(N+1)
Ab(1:N,1:N)=A
Ab(:,N+1)=b

do k=1,N-1
    elmax=dabs(Ab(k,k))
    id_max=k
  
 do i=k+1,n
      if (dabs(Ab(i,k))>elmax) then
         elmax=Ab(i,k)
         id_max=i
      end if          
    end do 
    vtemp1=Ab(k,:)
    vtemp2=Ab(id_max,:)   
    Ab(k,:)=vtemp2
    Ab(id_max,:)=vtemp1   

     do i=k+1,N  
     temp=Ab(i,k)/Ab(k,k)     
     Ab(i,:)=Ab(i,:)-temp*Ab(k,:)   
   end do
end do
Aup(:,:)=Ab(1:N,1:N)
bup(:)=Ab(:,N+1)
!调用用上三角方程组的回带方法
call uptri(Aup,bup,x,n)
end subroutine elgauss
subroutine uptri(A,b,x,N)
implicit real*8(a-z)
integer::i,j,N
real*8::A(N,N),b(N),x(N)
x(N)=b(N)/A(N,N)
!回带部分
do i=n-1,1,-1   
    x(i)=b(i)
   do j=i+1,N
    x(i)=x(i)-a(i,j)*x(j)
   end do
    x(i)=x(i)/A(i,i)
end do
end subroutine uptri
end module inv_mat
module  m_bfs
use inv_mat
contains
subroutine solve(x0,N,tol)
implicit real*8(a-z)
integer::i,n,itmax=50
real*8::x1(n),x0(n),y(n),f0(n),f1(n),dx(n)
real*8::v1(n),v2(n),v3(n)
real*8::H0(n,n),H1(n,n),df(n,n)
real*8::m1(n,n),m2(n,n),m3(n,n)
!itmax 最大允许迭代次数
!tol 误差容限
call jac(df,x0)
call inv(df,H0,N)
write(11,101)
write(12,102)
do i=1,itmax    
    !计算函数值
    call func(f0,x0)    
    !更新方程的根
    x1=x0-matmul(H0,f0)   
    dx=x1-x0   
    call func(f1,x1)    
    y=f1-f0    
    !这段程序求得miu
    v1=matmul(H0,y)
    t1=vdot(y,v1,n)
    t2=vdot(dx,y,N)
    u=1d0+t1/t2    
    !分子第一项
    m1=vvmat(dx,dx,n)*u    
    !分子第二项
    m2=vvmat(dx,y,n)
    m2=matmul(m2,H0)    
    !分子第三项
    v3=matmul(H0,y)
    m3=vvmat(v3,dx,n)    
    !分母
    t3=vdot(dx,y,n)    
    H1=(m1-m2-m3)/t3    
    H1=H0+H1    
    !把计算中的H矩阵写入文件
    write(12,103)i,H0    
    x0=x1
    H0=H1    
    !把计算迭代序列写入文件
    write(11,104)i,x0     
    !判断计算精度，当满足误差容限时 退出循环。
    dx2=dsqrt(dx(1)**2+dx(2)**2) 
 if (dx2<tol) exit 
end do
101 format(/,T16,'BFS方法计算序列',/)
102 format(/,T5,'BFS方法H矩阵序列为：',/)
103 format(T5,'iter=',I4,/,<N>(<N>F16.10/),/)
104 format(I4,3F16.10)
end subroutine solve
function vdot(a,b,N)
integer::i,N
real*8::a(n),b(n),vdot
vdot=0
do i=1,N  
  vdot=vdot+a(i)*b(i)
end do
end function vdot
function vvmat(a,b,n)
integer::n,i,j
real*8::a(n),b(n),vvmat(n,n)
do i=1,n
  do j=1,n
    vvmat(i,j)=a(i)*b(j)
  end do
end do
end function vvmat
subroutine func(f,x)
implicit real*8(a-z)
real*8::x(4),f(4),pi=3.141592653589793
real*8::xa=1.0,ya=1.0,za=1.0,xb=3.0,yb=1.0,zb=1.0
real*8::function_a=1.0,function_b=1.0,function_c=1.0,s1=2.0,s2=2.0
real*8::b11a=1.0, b12a=1.0, b13a=1.0, b21a=1.0, b22a=1.0, b23a=1.0, b31a=1.0, b32a=1.0, b33a=1.0
real*8::b11b=1.0, b12b=1.0, b13b=1.0, b21b=1.0, b22b=1.0, b23b=1.0, b31b=1.0, b32b=1.0, b33b=1.0
real*8::A(3), B(3)
A(1)=abs((b11a*(x(1)-xa)+b12a*(x(2)-ya)+b13a*(x(3)-za))/function_a)
A(2)=abs((b21a*(x(1)-xa)+b22a*(x(2)-ya)+b23a*(x(3)-za))/function_b)
A(3)=abs((b31a*(x(1)-xa)+b32a*(x(2)-ya)+b33a*(x(3)-za))/function_c)
B(1)=abs((b11b*(x(1)-xb)+b12b*(x(2)-yb)+b13b*(x(3)-zb))/function_a)
B(2)=abs((b21b*(x(1)-xb)+b22b*(x(2)-yb)+b23b*(x(3)-zb))/function_b)
B(3)=abs((b31b*(x(1)-xb)+b32b*(x(2)-yb)+b33b*(x(3)-zb))/function_c)
f(1)= 
f(2)= 
f(3)=
f(4)=
end subroutine func
subroutine  jac(df,x)
implicit real*8(a-z)
real*8::x(4),df(4,4)
real*8::xa=1.0,ya=1.0,za=1.0,xb=3.0,yb=1.0,zb=1.0
real*8::function_a=1.0,function_b=1.0,function_c=1.0,s1=2.0,s2=2.0
real*8::b11a=1.0, b12a=1.0, b13a=1.0, b21a=1.0, b22a=1.0, b23a=1.0, b31a=1.0, b32a=1.0, b33a=1.0
real*8::b11b=1.0, b12b=1.0, b13b=1.0, b21b=1.0, b22b=1.0, b23b=1.0, b31b=1.0, b32b=1.0, b33b=1.0
real*8::A(3), B(3)
A(1)=abs((b11a*(x(1)-xa)+b12a*(x(2)-ya)+b13a*(x(3)-za))/function_a)
A(2)=abs((b21a*(x(1)-xa)+b22a*(x(2)-ya)+b23a*(x(3)-za))/function_b)
A(3)=abs((b31a*(x(1)-xa)+b32a*(x(2)-ya)+b33a*(x(3)-za))/function_c)
B(1)=abs((b11b*(x(1)-xb)+b12b*(x(2)-yb)+b13b*(x(3)-zb))/function_a)
B(2)=abs((b21b*(x(1)-xb)+b22b*(x(2)-yb)+b23b*(x(3)-zb))/function_b)
B(3)=abs((b31b*(x(1)-xb)+b32b*(x(2)-yb)+b33b*(x(3)-zb))/function_c)
df(1,1)=
df(1,2)=
df(1,3)=
df(1,4)=
df(2,1)=
df(2,2)=
df(2,3)=
df(2,4)=
df(3,1)=(s1-s2)*(A(1)**s2+A(2)**s2)**(s1/s2-2.0d0)*(b13a/function_a*A(1)**(s2-1.0d0)+b23a/function_b*A(2)**(s2-1.0d0))*(b11a &
    &  /function_a*A(1)**(s2-1.0d0)+b21a/function_b*A(2)**(s2-1.0d0)) &
    & +(A(1)**s2+A(2)**s2)**(s1/s2-1.0d0)*(b11a*b13a/function_a**2.0d0*(s2-1.0d0)*A(1)**(s2-2.0d0)+b21a*b23a/fuction_b**2.0d0 &
    & *(s2-1.0d0)*A(2)**(s2-2.0d0))+b31a*b33a/function_c**2.0d0*(s1-1.0d0)*A(3)**(s1-2.0d0)+x(4)*((s1-s2)*(B(1)**s2+B(2)**s2) &
    & **(s1/s2-2.0d0)*(b13b/function_a*B(1)**(s2-1.0d0)+b23b/function_b*B(2)**(s2-1.0d0))*(b11b &
    &  /function_a*B(1)**(s2-1.0d0)+b21b/function_b*B(2)**(s2-1.0d0)) &
    & +(B(1)**s2+B(2)**s2)**(s1/s2-1.0d0)*(b11b*b13b/function_a**2.0d0*(s2-1.0d0)*B(1)**(s2-2.0d0)+b21b*b23b/fuction_b**2.0d0 &
    & *(s2-1.0d0)*B(2)**(s2-2.0d0))+b31b*b33b/function_c**2.0d0*(s1-1.0d0)*B(3)**(s1-2.0d0))
df(3,2)=(s1-s2)*(A(1)**s2+A(2)**s2)**(s1/s2-2.0d0)*(b12a/function_a*A(1)**(s2-1.0d0)+b22a/function_b*A(2)**(s2-1.0d0))*(b13a &
    &  /function_a*A(1)**(s2-1.0d0)+b23a/function_b*A(2)**(s2-1.0d0)) &
    & +(A(1)**s2+A(2)**s2)**(s1/s2-1.0d0)*(b13a*b12a/function_a**2.0d0*(s2-1.0d0)*A(1)**(s2-2.0d0)+b23a*b22a/fuction_b**2.0d0 &
    & *(s2-1.0d0)*A(2)**(s2-2.0d0))+b33a*b32a/function_c**2.0d0*(s1-1.0d0)*A(3)**(s1-2.0d0)+x(4)*((s1-s2)*(B(1)**s2+B(2)**s2) &
    & **(s1/s2-2.0d0)*(b12b/function_a*B(1)**(s2-1.0d0)+b22b/function_b*B(2)**(s2-1.0d0))*(b13b &
    &  /function_a*B(1)**(s2-1.0d0)+b23b/function_b*B(2)**(s2-1.0d0)) &
    & +(B(1)**s2+B(2)**s2)**(s1/s2-1.0d0)*(b13b*b12b/function_a**2.0d0*(s2-1.0d0)*B(1)**(s2-2.0d0)+b23b*b22b/fuction_b**2.0d0 &
    & *(s2-1.0d0)*B(2)**(s2-2.0d0))+b33b*b32b/function_c**2.0d0*(s1-1.0d0)*B(3)**(s1-2.0d0))
df(3,3)=(s1-s2)*(A(1)**s2+A(2)**s2)**(s1/s2-2.0d0)*(b13a/function_a*A(1)**(s2-1.0d0)+b23a/function_b*A(2)**(s2-1.0d0))**2.0d0 &
    & +(A(1)**s2+A(2)**s2)**(s1/s2-1.0d0)*((b13a/function_a)**2.0d0*(s2-1.0d0)*A(1)**(s2-2.0d0)+(b23a/fuction_b)**2.0d0 &
    & *(s2-1.0d0)*A(2)**(s2-2.0d0))+(b33a/function_c)**2.0d0*(s1-1.0d0)*A(3)**(s1-2.0d0)+x(4)*((s1-s2)*(B(1)**s2+B(2)**s2) &
    & **(s1/s2-2.0d0)*(b13b/function_a*B(1)**(s2-1.0d0)+b23b/function_b*B(2)**(s2-1.0d0))**2.0d0+(B(1)**s2+B(2)**s2)**(s1/s2 &
    & -1.0d0)*((b13b/function_a)**2.0d0*(s2-1.0d0)*B(1)**(s2-2.0d0)+(b23b/fuction_b)**2.0d0 *(s2-1.0d0)*B(2)**(s2-2.0d0)) &
    & +(b33b/function_c)**2.0d0*(s1-1.0d0)*B(3)**(s1-2.0d0))
df(3,4)=(B(1)**s2+B(2)**s2)**(s1/s2-2.0d0)*(b13b/function_a*B(1)**(s2-1.0d0)+b23b/function_b*B(2)**(s2-1.0d0)) &
    & +b33b/function_c*B(3)**(s1-1.0d0)
df(4,1)=(B(1)**s2+B(2)**s2)**(s1/s2-2.0d0)*(b11b/function_a*B(1)**(s2-1.0d0)+b21b/function_b*B(2)**(s2-1.0d0)) &
    & +b31b/function_c*B(3)**(s1-1.0d0)
df(4,2)=(B(1)**s2+B(2)**s2)**(s1/s2-2.0d0)*(b12b/function_a*B(1)**(s2-1.0d0)+b22b/function_b*B(2)**(s2-1.0d0)) &
    & +b32b/function_c*B(3)**(s1-1.0d0)
df(4,3)=(B(1)**s2+B(2)**s2)**(s1/s2-2.0d0)*(b13b/function_a*B(1)**(s2-1.0d0)+b23b/function_b*B(2)**(s2-1.0d0)) &
    & +b33b/function_c*B(3)**(s1-1.0d0)
df(4,4)=0.0d0
end subroutine jac
end module m_bfs
program main
use  inv_mat
use  m_bfs
implicit real*8(a-z)
integer::N=4
real*8::x0(4)
open(unit=11,file='result.txt')
open(unit=12,file='Hmatrix.txt')
x0=(/2.0d0,1.0d0,1.0d0,1.0d0/)
!1d-8  表示允许的误差容限
call solve(x0,n,1d-8)
end program main

想问下大家，这个程序有什么错误吗？为什么最后的结果总是NAN呢？

pasuka

矩阵线性运算为啥不用OpenBLAS或者MKL呢？
自编代码，尤其是徐*良的代码坑人不浅。。。。

无奈的小强

pasuka 发表于 2016-8-24 13:26
矩阵线性运算为啥不用OpenBLAS或者MKL呢？
自编代码，尤其是徐*良的代码坑人不浅。。。。 ...

能问下您，您有比较成熟的模板可以求解我这种复杂非线性方程组的吗？

pasuka

无奈的小强 发表于 2016-8-24 15:09
能问下您，您有比较成熟的模板可以求解我这种复杂非线性方程组的吗？

1、非线性方程求解目前没有理想的通用方法
2、强非线性问题可以试试廖世俊老师的同伦方法
3、科研性质的原理探索与验证，还是MATLAB事半功倍

无奈的小强

pasuka 发表于 2016-8-24 15:18
1、非线性方程求解目前没有理想的通用方法
2、强非线性问题可以试试廖世俊老师的同伦方法
3、科研性质的 ...

好的，谢谢您~！

kerb

其实这个问题不难，只是你加了一个限制条件，“拟牛顿BFS方法”，不知道你对这个方法是否熟悉，如果不熟悉，你很难找到哪里出错
你的问题实际上是求解一个非线性方程组：有4个未知数，4个非线性方程组成方程组，首先你最好在纸上手工写写画画，搞明白再动手写代码F(X)=(f_1(X),f_2(X),f_3(X),f_4(x))^T;
其中X=(x_1,x_2,x_3,x_4)
第一步你要构造一个目标函数：OBJFUNC(X)=\frac{1}{2}\sum_{i=1}^4 f_i^2(X)，如果方程组F(X)=0有解，对应OBJFUNC(X)的最小值
第二步你要把F(X)taylor展开：F(X_0+\Delta X)=F(X_0)+\nabla F \Delta X + \cdots，其中\nabla F叫雅克比矩阵，一般用J表示，对你的这个问题是一个4X4的矩阵
第三步，求解OBJFUNC(X)的最小值，这个时候，其实4个方程4个变量，你应该可以解析计算出HESSIAN矩阵的，考虑到现在的孩子们大部分时间用来看手机，脑筋不太够用，不会算也是很有可能的，那么一个讨巧的办法是可以用的，就是利用第二步的那个泰勒展开：同时假设F(X_0+\Delta X)，这样就有F(X_0)+\nabla F \Delta X=0或者F+J \Delta X=0，如果\Delta X不大，可以证明海赛矩阵与J^TJ几本差不多，这样问题就可以转化为求解：J^TJ \Delta X\approx-J^T F
第四步你就可以跌代求解了：需要明确的是这里的若干符号的具体含义，实际上J\rightarrow J(X_n),F\rightarrow F(x_n), \Delta X=X_{n+1}-X_n
J^TJ是对称的，但不一定是正定的，很多时候是奇异或接近奇异，跌代过程往往没法进行下去，Marquardt-Levenber搞了一个权宜之计，就是每轮跌代过程中，给J^TJ的对角元素加上一个很小的正数\mu，防止J^TJ奇异或者条件数很大，每轮跌代过程中都需要修改\mu，以防止\mu的阻碍作用(这种方法实际上还有一个名字：阻尼牛顿法或阻尼最小二乘法)，如何修改\mu也很有学问，一个简单偷懒的办法可以这样：根据每轮OBJFUNC(X)的大小来调整\mu，具体就是本轮迭代中如果OBJFUNC(X)比上一轮的大，\mu=\mu*4.0，反之\mu=\mu/2.0跌代最初的时候令\mu=1.E-3。
如果OBJFUNC(X)非常小了，也就是满足你的精度要求了，就可以终止跌代

kerb

拟牛顿BFS方法是另一种防止跌代过程中海赛矩阵奇异或者条件数很大的办法，里面有很多矩阵或者矢量相乘的运算，需要细心，搞不好就出错

无奈的小强

kerb 发表于 2016-8-30 14:29
其实这个问题不难，只是你加了一个限制条件，“拟牛顿BFS方法”，不知道你对这个方法是否熟悉，如果不熟悉 ...

好的，谢谢您，我再试一下。