Fortran学习19：数值计算方法2

最新推荐文章于 2022-10-17 11:07:04 发布
最新推荐文章于 2022-10-17 11:07:04 发布
阅读量2.5k 收藏 13
点赞数 1
分类专栏： Fortran 文章标签： 学习 开发语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/faltas/article/details/127315946
版权

老感觉昨天自己最后写的那个python牛顿法求解，老觉得不太对，，其实后面我想学完Fortran之后，再跟着matlab书把数据处理再学一下，因为Fortran那本书讲的东西实在是太少了。，而且听师兄说，我们不用Fortran绘图，后面还是matlab和python来绘图。

今天先不说了，直接开始矩阵运算，其实矩阵的加减法没啥好说的，就是相同位置进行加减运算就行了，主要是乘法，Fortran库里有一个matmul()函数来计算，但是今天还是先自己敲吧。

program main
    implicit none
    integer::a(3,3)
    integer::b(3,3)
    integer::c(3,3)
    integer i,j,k
    data a / 1,2,3,4,5,6,7,8,9 /
    data b /10,11,12,13,14,15,16,17,18/
    ! 加法
    ！可以直接使用c=a+b，下面是为了展示运算过程
    do i=1,3
        do j=1,3
            c(i,j)=a(i,j)+b(i,j)
        end do
    end do
    ! 减法
    ！可以直接使用c=a-b，下面是为了展示运算过程
    do i=1,3
        do j=1,3
            c(i,j)=a(i,j)-b(i,j)
        end do
    end do
    ! 乘法
    do i=1,3
        do j=1,3
            c(i,j)=0
            do k=1,3
                c(i,j)=c(i,j)+a(i,k)*b(k,j)
            end do 
        end do
    end do
    ! 再试试Fortran自带的函数matmul
    c=MATMUL(a,b)
    write(*,"(I5,I5,I5/I5,I5,I5/I5,I5,I5)")c
    stop "over"
end program

放心，没问题的，ctrl+c，ctrl+v之后完全可以运行。

其实也可以在主程序外面把这几个运算方法写成function或者subroutine，我觉得都还行，就是个传参嘛。反正我是没有用书上的代码，习惯了，不喜欢整个程序纯粹用大写字母。

1、三角矩阵：上三角和下三角，上三角对角线以下的数全部是0，相应的，下三角对角线以上的全部是零：

跟着书上的代码敲一遍吧，我线代的储备主要还是考研时候听网课攒的点底子(大二学线代的时候才考了85+，有点低)

module LinearAlgebra
    implicit none
    contains
        subroutine output(matrix)
            implicit none
            integer::m,n
            real::matrix(:,:)
            integer::i
            CHARACTER(len=20)::for='(??(1x,f6.3))'
            m=size(matrix, dim=1)
            n=size(matrix, dim=2)
            WRITE(for(2:3),'(I2)')N
            do i=1,N
                write(*,Fmt=for)matrix(i,:)
            end do
            return 
        end subroutine 
        ! 求解上三角矩阵：
        subroutine upper(matrix)
            implicit none
            integer::m,n
            real::matrix(:,:)
            integer::i,j
            real::e
            m=size(matrix, dim=1)
            n=size(matrix, dim=2)
            do I=1,n-1
                do j=i+1,m
                    e=matrix(j,i)/matrix(i,i)
                    matrix(j,i:m)=matrix(j,i:m)-matrix(i,i:m)*e
                end do
            end do
            return 
        end subroutine
        ! 求解下三角矩阵：
        subroutine lower(matrix)
            implicit none
            integer::m,n
            real::matrix(:,:)
            integer::i,j
            real::e
            m=size(matrix, dim=1)
            n=size(matrix, dim=2)
            do I=n,2,-1
                do j=i-1,1,-1
                    e=matrix(j,i)/matrix(i,i)
                    matrix(j,1:i)=matrix(j,1:i)-matrix(i,1:i)*e
                end do
            end do
            return 
        end subroutine
end module
program main
    use LinearAlgebra
    implicit none
    INTEGER,PARAMETER::n=3
    real::a(n,n),b(n,n)
    data a / 1,2,1,3,2,3,2,3,4 /
    write(*,*)"Matrix A:"
    CALL output(a)
    b=a
    write(*,*)"Upper:"
    CALL upper(b)
    CALL output(b)
    b=a
    write(*,*)"lower:"
    CALL upper(b)
    CALL output(b)
    stop "over"
end program

输出：

 学习上下三角矩阵的目的就是，可以应用在determinant、gauss-jordan发求解联立式、求逆矩阵等。

下午如果有时间的话我就学着写写determinant、gauss-jordan和求逆矩阵如果没时间的话那就算了。

哈哈，下午有时间！

Determinant：就是把矩阵转化成上三角矩阵，然后求行列式值，把这段代码加到module里面就可以调用。

        real function Determinant(matrix)
            real::matrix(:,:)
            real,ALLOCATABLE::ma(:,:)
            integer::i,n
            n=size(matrix, dim=1)
            ALLOCATE(ma(n,n))
            ma=matrix
            CALL upper(ma)
            Determinant=1.0
            do i=1,n
                Determinant=Determinant*ma(i,i)
            end do
        end function

 Gauss-Jordan法求联立方程式：

        subroutine Gauss_Jordan(a,s,ans)
            real::a(:,:)
            real::s(:)
            real::ans(:)
            real,ALLOCATABLE::b(:,:)
            integer :: i
            integer :: n 
            n=size(a, dim=1)
            ALLOCATE(b(n,n))
            b=a
            ans=s
            call upper1(b,ans,n)
            call lower1(b,ans,n)
            do i=1,n
                ans(i)=ans(i)/b(i,i)
            end do
            RETURN
        end subroutine
        subroutine output1(m,s)
            implicit none
            real::m(:,:),s(:)
            integer i,j,n
            n=size(m, dim=1)
            do i=1,n
                write(*,"(1x,f5.2,a1)",advance="no")m(i,1),'a'
                do j=2,n
                    if(m(i,j)<0)then
                        write(*,"('-',f5.2,a1)",advance="no")-m(i,j),char(64+j)
                    else
                        write(*,"('+',f5.2,a1)",advance="no")m(i,j),char(64+j)
                    end if
                end do 
            end do 
        end subroutine
        subroutine upper1(matrix,s,n)
            implicit none
            integer::n
            real::matrix(n,n)
            real::s(n)
            integer::i,j
            real::e
            do I=1,n-1
                do j=i+1,n
                    e=matrix(j,i)/matrix(i,i)
                    matrix(j,i:n)=matrix(j,i:n)-matrix(i,i:n)*e
                    s(j)=s(j)-s(i)*e
                end do
            end do
            return 
        end subroutine
        ! 求解下三角矩阵：
        subroutine lower1(matrix,s,n)
            implicit none
            integer::n
            real::matrix(:,:)
            real::s(n)
            integer::i,j
            real::e
            do I=n,2,-1
                do j=i-1,1,-1
                    e=matrix(j,i)/matrix(i,i)
                    matrix(j,1:n)=matrix(j,1:n)-matrix(i,1:n)*e
                    s(j)=s(j)-s(i)*e
                end do
                write(*,"('=',f8.4)")s(i)
            end do
            return 
        end subroutine

上面两段代码加到module里面就好，直接运行。

 Gauss_Jordan是用来求解方程组的。

逆矩阵，实际上书上介绍的方法就是在要求逆矩阵的矩阵后面加上一个单位矩阵，然后把前者处理成单位矩阵，前者每一次处理的时候，后者要跟着进行相同的操作，之后得到的就是逆矩阵。

        subroutine inverse(a,b)
            implicit none
            real::a(:,:),b(:,:)
            real,ALLOCATABLE::c(:,:)
            integer i,j,n
            n=size(a, dim=1)
            allocate(c(n,n))
            forall(i=1:n,j=1:n,i==j)b(i,j)=1.0
            forall(i=1:n,j=1:n,i/=j)b(i,j)=0.0
            c=a
            call upper1(c,b,n)
            call lower1(c,b,n)
            forall(i=1:n)b(i,:)=b(i,:)/b(i,i)
            return
        end subroutine

 

这下求出来的就是matrix的逆矩阵啦，同样的，上面的subroutine直接加到module里面就可以调用它了。

至于使用嘛，就看你到底要处理什么问题了，比如说求解下面的对角矩阵啊，解方程组之类的，都可以，具体情况具体看待。

\begin{bmatrix} 2 3 0 0 0 0 0\\ 1 2 3 0 0 0 0\\ 0 1 2 3 0 0 0 \\ 0 0 1 2 3 0 0 \\ 0 0 0 1 2 3 0\\ 0 0 0 0 1 2 3\\ 0 0 0 0 0 1 2 \end{bmatrix}

ok，把今天的完整代码附上：

module LinearAlgebra
    implicit none
    contains
        !  Determinant矩阵
        real function Determinant(matrix)
            real::matrix(:,:)
            real,ALLOCATABLE::ma(:,:)
            integer::i,n
            n=size(matrix, dim=1)
            ALLOCATE(ma(n,n))
            ma=matrix
            CALL upper(ma)
            Determinant=1.0
            do i=1,n
                Determinant=Determinant*ma(i,i)
            end do
        end function 
        
        ! 求解联立方程组
        subroutine Gauss_Jordan(a,s,ans)
            real::a(:,:)
            real::s(:)
            real::ans(:)
            real,ALLOCATABLE::b(:,:)
            integer :: i
            integer :: n 
            n=size(a, dim=1)
            ALLOCATE(b(n,n))
            b=a
            ans=s
            call upper1(b,ans,n)
            call lower1(b,ans,n)
            do i=1,n
                ans(i)=ans(i)/b(i,i)
            end do
            RETURN
        end subroutine

        ! 逆矩阵求解
        subroutine inverse(a,b)
            implicit none
            real::a(:,:),b(:,:)
            real,ALLOCATABLE::c(:,:)
            integer i,j,n
            n=size(a, dim=1)
            allocate(c(n,n))
            forall(i=1:n,j=1:n,i==j)b(i,j)=1.0
            forall(i=1:n,j=1:n,i/=j)b(i,j)=0.0
            c=a
            call upper1(c,b,n)
            call lower1(c,b,n)
            forall(i=1:n)b(i,:)=b(i,:)/b(i,i)
            return
        end subroutine

        ! 输出
        subroutine output(matrix)
            implicit none
            integer::m,n
            real::matrix(:,:)
            integer::i
            CHARACTER(len=20)::for='(??(1x,f6.3))'
            m=size(matrix, dim=1)
            n=size(matrix, dim=2)
            WRITE(for(2:3),'(I2)')N
            do i=1,N
                write(*,Fmt=for)matrix(i,:)
            end do
            return 
        end subroutine 

        ! 求解上三角矩阵：
        subroutine upper(matrix)
            implicit none
            integer::m,n
            real::matrix(:,:)
            integer::i,j
            real::e
            m=size(matrix, dim=1)
            n=size(matrix, dim=2)
            do I=1,n-1
                do j=i+1,m
                    e=matrix(j,i)/matrix(i,i)
                    matrix(j,i:m)=matrix(j,i:m)-matrix(i,i:m)*e
                end do
            end do
            return 
        end subroutine
        
        ! 求解下三角矩阵：
        subroutine lower(matrix)
            implicit none
            integer::m,n
            real::matrix(:,:)
            integer::i,j
            real::e
            m=size(matrix, dim=1)
            n=size(matrix, dim=2)
            do I=n,2,-1
                do j=i-1,1,-1
                    e=matrix(j,i)/matrix(i,i)
                    matrix(j,1:i)=matrix(j,1:i)-matrix(i,1:i)*e
                end do
            end do
            return 
        end subroutine

        ! 输出1：供gauss——jordan使用
        subroutine output1(m,s)
            implicit none
            real::m(:,:),s(:)
            integer i,j,n
            n=size(m, dim=1)
            do i=1,n
                write(*,"(1x,f5.2,a1)",advance="no")m(i,1),'a'
                do j=2,n
                    if(m(i,j)<0)then
                        write(*,"('-',f5.2,a1)",advance="no")-m(i,j),char(64+j)
                    else
                        write(*,"('+',f5.2,a1)",advance="no")m(i,j),char(64+j)
                    end if
                end do 
            end do 
        end subroutine

        ! 求解上三角矩阵1
        subroutine upper1(matrix,s,n)
            implicit none
            integer::n
            real::matrix(n,n)
            real::s(n)
            integer::i,j
            real::e
            do I=1,n-1
                do j=i+1,n
                    e=matrix(j,i)/matrix(i,i)
                    matrix(j,i:n)=matrix(j,i:n)-matrix(i,i:n)*e
                    s(j)=s(j)-s(i)*e
                end do
            end do
            return 
        end subroutine
        
        ! 求解下三角矩阵1：
        subroutine lower1(matrix,s,n)
            implicit none
            integer::n
            real::matrix(:,:)
            real::s(n)
            integer::i,j
            real::e
            do I=n,2,-1
                do j=i-1,1,-1
                    e=matrix(j,i)/matrix(i,i)
                    matrix(j,1:n)=matrix(j,1:n)-matrix(i,1:n)*e
                    s(j)=s(j)-s(i)*e
                end do
                write(*,"('=',f8.4)")s(i)
            end do
            return 
        end subroutine
end module
program main
    use LinearAlgebra
    implicit none
    INTEGER,PARAMETER::n=3
    real::a(n,n),b(n,n)
    real::a1(n,n),s(n),ans(n)
    integer i
    data a / 1,2,1,3,2,3,2,3,4 /
    data a1 / 1,2,3,4,5,6,7,8,8 /
    data s / 12,15,17 /

    write(*,*)"equation"
    CALl output1(a,s)
    call Gauss_Jordan(a,s,ans)
    write(*,*)"ans="
    
    do i=1,n
        write(*,"(1x,a1,'=',f8.4)")char(64+i),ans(i)
    end do
    
    write(*,*)"matrix:"
    call output(a)
    call inverse(a,b)
    
    write(*,*)"inverse"
    call output(b)
    
    write(*,*)"Matrix A:"
    CALL output(a)

    b=a
    write(*,*)"Upper:"
    CALL upper(b)
    CALL output(b)
    
    b=a
    write(*,*)"lower:"
    CALL upper(b)
    CALL output(b)
    write(*,"('det(A)=',F6.2)")Determinant(a)
    stop "over"
end program

计算技术【40】
基础数学研究
09-26 869
FORTRAN是英文“FORmulaTRANslator”的缩写，译为“公式翻译器”，它是世界上最早出现的计算机高级程序设计语言，广泛应用于科学和工程计算领域。FORTRAN语言以其特有的功能在数值、科学和工程计算领域发挥着重要作用。Fortran语言的最大特性是接近数学公式的自然描述，在计算机里具有很高的执行效率。易学，语法严谨。可以直接对矩阵和复数进行运算，这一点类似matlab。自诞生以来广泛地应用于数值计算领域，积累了大量高效而可靠的源程序。
fortran计算实例
05-30
fortran计算实例
Fortran学习18：数值计算方法1
faltas的博客
10-13 2043
Fortran学习：数值分析（二分法、割线法、牛顿法）
Fortran学习22：数值计算方法5
faltas的博客
10-17 882
Fortran学习
Fortran矩阵运算
热门推荐
heyzol的博客
11-09 3万+
Fortran中自带了一个进行矩阵乘法运算的函数 MATMUL，用来计算矩阵乘法。 矩阵在fortran中，一个3*3的数组a(3,3)所表示的矩阵如下 program main implicit none real a(3, 3), b(3, 3), c(3, 3) integer m, n, k integer i, j a = reshap
数值分析—全主元消去法解线性方程组—FORTRAN程序
qq_43630847的博客
11-12 1793
数值分析—全主元消去法解线性方程组—FORTRAN程序 program main implicit none real8,dimension(:
Fortran95学习笔记：从C语言到科学计算
笔记作者强调了Fortran在数值计算和科学数据处理中的重要性，并对比了Fortran77和Fortran90的主要差异。 Fortran，全称Formula Translator或Translation，是一种早期的高级编程语言，特别设计用于数值计算和科学...
基础Fortran代码包：完全可运行的编程学习资源
Fortran语言在数值计算方面有广泛的应用，例如在解决线性代数问题、微分方程求解、傅里叶变换等领域。它提供了高级的数学函数库，支持复杂数学运算的实现。 7. 程序的调试和优化 由于Fortran程序通常用于科学计算，...
Fortran95学习笔记：彭国伦《Fortran95程序设计》要点解析
自IBM在早期开发以来，它就以其强大的数组操作和对数值计算的支持，在科研和工程领域中广泛应用。Fortran的发展历程中产生了多个版本，其中Fortran77和Fortran90是使用较为广泛的两个。 Fortran77是早期的版本，它...
Fortran95学习笔记：彭国伦著书要点解析
在深入探讨之前，首先理解Fortran是一种古老而强大的编程语言，最初设计用于科学计算和数值分析，其名称是Formula Translator的缩写。Fortran的最新版本，如Fortran90，引入了许多改进，使得它更现代、更易于编写和...
Fortran PowerStation 4.0：科学计算与工程领域的高效编译器
FORTRAN语言的特点是结构化编程，支持数组运算和复数运算，非常适合用于大规模数值计算。 2. Fortran PowerStation 4.0:Fortran PowerStation 4.0是一款针对FORTRAN语言开发的编译器，由Microsoft公司发布。它具备...
fortran数值分析算法
09-11
数值分析的fortran算法大全
fortran标准函数库
11-21
fortran标准函数库，包含字符运算、矩阵计算等fortran编译器自带算法。
fortran 矩阵运算
05-16
fortran 矩阵运算，进行矩阵的加、乘、求逆、转置，FORTRAN源码，经过验证
Fortran计算矩阵2范数
chd_lkl的博客
03-19 4299
Intel fortran编译器自带的MKL函数库中，只有计算矩阵1范数和无穷范数的程序，所以本文给出fortran使用MKL函数库计算矩阵2范数的示例。 而矩阵的2范数的值近似于max(svd(A)) 运行环境：win10+vs2019+ivf2020 代码如下： program norm2 use lapack95 implicit none ...
Fortran入门教程(四)——数据运算
Sumbrella_的博客
12-10 6159
#数据运算 1. 运算符 fortran中最基本的五种运算符为 + 加 - 减 * 乘 / 除 ** 乘方 mod(a, b) a对b取余 优先级: 乘法 > 乘除 > 加减, 同一级别从左到右依次运算。 2. 运算关系符 记号表达式 字母表达式 含义 < .lt. 小于 <= .le. 小于等于 == .eq. 等于 /= .ne. 不等于 > .gt. 大于 >= .ge. 大于等于 在使用时，可以使用记号表达式也可
fortran：关于矩阵运算的matmul计算速度过慢解决方案
chd_lkl的博客
06-02 5491
下面给出两个方阵相乘的示例代码. Z(m,L) = X(m,n) * Y(n,L). 其中X,Y,Z也可以是长阵 program main use omp_lib !// 使用omp_lib并行 use ifport, only : fdate !// 使用fdate进行计算时间统计 implicit none integer,...
Fortran几个函数（DOT_PRODUCT+MATMUL+TRANSPOSE+RESHAPE）
Jiwei的专栏
01-13 1万+
DOT_PRODUCT函数（点积）     DOT_PRODUCT(VECTOR_A, VECTOR_B) computes the dot product multiplication of two vectors VECTOR_A and VECTOR_B. The two vectors may be either numeric or logical and must b
矩阵每一行重复_Fortran如何实现矩阵与向量的乘法运算
weixin_31578739的博客
12-24 1079
矩阵是二维数组，而向量是一维数组，内置函数matmul不能实现矩阵与向量的乘法运算。在这一点Fortran不如matlab灵活。Fortran如何实现矩阵与向量的乘法运算，现有以下三种方法供参考。一)将一维数组看作二维数组的退化形式，比如a(3)可以看作a(3,1)或者a(1,3)，这样就可以用matmul函数计算了。二)用spread函数将一维数组扩展成二维数组，同样可用matmul函...
评论 8
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值