Gram-Schmidt计算正交基Fortran代码

0. 引言

  Gram-Schmidt过程是一种用于将一组线性无关的向量转换为一组正交或标准正交向量的算法。这个过程通常用于线性代数和数值计算中。

  假设我们有一组线性无关的向量 $v_1,v_2,\dots ,v_n$ 。Gram-Schmidt过程的目标是生成一组正交向量 $q_1,q_2,\dots ,q_n$，它们仍然在原始向量的线性空间中，但是彼此之间是正交的（或标准正交的，即彼此之间的内积为0，且每个向量的范数为1）。

Gram-Schmidt过程的步骤：

a) 初始化：

  将第一个向量作为新的正交向量的第一个成员：
$$q_1=v_1$$

b 递归构建正交向量：
  对于每个 $i=2,3,\dots ,n$，执行以下步骤：

• 正交化处理：

  计算当前向量 $v_i$ 在前 $i-1$ 个已经确定的正交向量 $[q_1,q_2,\dots ,q_{i-1}]$ 上的投影。投影的计算方式为：
$$\text{proj}{q}_{i-1}{v}_i=\frac{⟨v_i,q_{i-1}⟩}{⟨q_{i-1},q_{i-1}⟩}{q}_{i-1}$$
  其中 $(⟨\cdot ,\cdot ⟩)$ 表示内积（或点积）。

• 更新向量：

  将投影向量从 (v_i) 中减去，得到一个新的向量 $(v_i^{\prime })$：
$$v_i^{\prime }=v_i-\text{proj}{q}_{i-1}(v_i)$$

• 正则化：

  将 $v_i^{\prime }$ 归一化为单位向量 $q_i$：

$$q_i=\frac{v_i^{\prime }}{|v_i^{\prime }|}$$
  其中 $|v_i^{\prime }|$表示向量 $v_i^{\prime }$ 的范数（或长度）。

• 输出结果：

  最终得到一组正交向量 $q_1,q_2,\dots ,q_n$，它们满足以下两个条件：
    (1) $q_i\cdot q_j=0$ 对于所有 $i\ne j$，即彼此正交。
    (2) $|q_i|=1$ 对于所有 $i$，即每个向量的范数为1。

  这样，Gram-Schmidt过程能够有效地将一个线性无关的向量集合转换为一组正交或标准正交的向量集合，为许多数学和计算问题提供了便利。

  本篇用fortran实现了Gram-Schmidt正交化过程，过程参考了matlab的相关资料。

1. Fortran实现Gram-Schmidt过程

  下面是示例程序，用简单案例测试Gram-Schmidt方法计算过程。

%%调用主程序
 program t
 use math
 real(8),allocatable :: aa(:,:),oaa(:,:)
 
 allocate(aa(3,3))
 aa(1,:) = [8d0,1d0,8d0];
 aa(2,:) = [8d0,1d0,4d0]   
 aa(3,:) = [6d0,0d0,5d0]   
    
 oaa = orths(aa) 
    
 print *,oaa
 
end program

执行结果：

%% 主调函数，输入向量 输出正交基向量
module math2
    implicit none

    contains

    function orths(A)result(ou)

    real(8),allocatable,intent(in) :: A(:,:)

    real(8),dimension(3) :: bt,at,at2,bt2
    integer :: N, i,ind
    real(8),allocatable,dimension(:,:) :: b,ou

    N = size(A,1)
    allocate( b(N,3) )
    do i = 1, N
        if(i == 1)then
            b(i,:) = A(i,:)
            cycle
        endif

        at(:) = A(i,:)
        at2 = at
        ind = 0
        do
            ind = ind + 1
            if(ind >= i)then
                exit
            endif
            bt(:) = b(ind,:)
            at2 = at2 - ( dot_product(at, bt)/dot_product(bt,bt) )* bt
        enddo
        b(i,:) = at2(:)
    enddo
    allocate(ou(N,3))
    ou = 0d0

    !> 归一化
    do i = 1, N
        bt2(:) = b(i,:)
        ou(i,:) = norm( bt2 )
    enddo

    contains

    ! norm
    function norm(a)result(b)

    real(8),intent(in) :: a(:)
    real(8) :: b(3)

    real(8) :: no

    no = sqrt( a(1)**2 + a(2)**2 + a(3)**2 )
    b = a / no
    end function

    end function

end module

2. 结语

   本篇分享了基于fortran的Gram-Schmidt方法，能够求解输入矩阵的正交基。希望有所帮助

😜
😜😜
😜😜😜😜