Fortran Coder

查看: 1051|回复: 4
打印 上一主题 下一主题

[非线性] 出现四个问题总是无法解决,求帮帮忙

[复制链接]

22

帖子

11

主题

0

精华

入门

F 币
97 元
贡献
54 点
跳转到指定楼层
楼主
发表于 2023-11-30 10:39:34 | 只看该作者 |只看大图 回帖奖励 |倒序浏览 |阅读模式
[Fortran] 纯文本查看 复制代码
Module m_gauss
Contains
Subroutine lineq(A,b,x,N)!!解线性方程组 高斯列主元消去法Ax=b
  Implicit real*8(a-z)
  Integer::i,k,N
  Integer::id_max !主元素符号
  real*8::A(N,N),b(N),x(N)!!A(N,N)系数矩阵,b(N)右向量,N方程维数,x方程的根
  real*8::Aup(N,N),bup(N)
  real*8::Ab(N,N+1)!Ab为增广矩阵[Ab]
  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
!!!查找主元素(不是赋值最大元素给elmax,而是为了找出最大元素对应的标号)
      do i=k+1,n
       if (dabs(Ab(i,k))>elmax) then
           elmax=Ab(i,k)
           id_max=i
       end if
      end do
!!!完成查找最大元素,查找完成以后与第k行交换,其他不变
   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 lineq

Subroutine uptri(A,b,x,N)!!A(N,N)系数矩阵,b(N)右向量,N方程维数,x方程的根
Implicit real*8(a-z)
Integer::i,j,k
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 m_gauss

Module m_newton
    Contains
    Subroutine solve()!!计算非线性方程组,需要用到线性方程组的解法,这里用主元消元法;需要准备函数文件和偏导数文件
    use m_gauss
    Implicit real*8(a-z)
    Integer::I,itmax=50 !itmax最大允许迭代次数
    Integer::N=2  !!N方程维数  
    real*8::x(2),f(2),dx(2)
    real*8::df(2,2)!!df偏导数矩阵
    open(unit=11,file='result.txt')
    write(11,101)
101 format(/,T6,'牛顿法计算非线性方程组迭代序列',/)    
    x=(/2d0,2d0/)
    tol=1d-8
    do i=1,itmax
        call func(f,x)
        call jac(df,x)
        call lineq(df,-f,dx,N)
        x=x+dx
        write(11,102)i,x
102  format(I5,2F16.10)       
     dx2=dsqrt(dx(1)**2+dx(2)**2)!判断计算精度,当满足误差容限时退出循环。deqrt是 平方根
     if(dx2<tol)exit
    end do
    end Subroutine solve

    Subroutine func(f,x)!f方程函数,x自变量.
    Implicit real*8(a-z)
    real*8::x(2),f(2)
    f(1)=6*x(1)**3+x(1)*x(2)-3*x(2)**3-4
    f(2)=x(1)**2-18*x(1)*x(2)**2+16*x(2)**3+1
    end Subroutine func

    Subroutine jac(df,x)
    Implicit real*8(a-z)
    real*8::x(2),df(2,2)!!偏导数矩阵
    df(1,1)=18*x(1)**2+x(2)
    df(2,1)=2*x(1)-18*x(2)**2
    df(1,2)=x(1)-9*x(2)**2
    df(2,2)=-36*x(1)*x(2)+48*x(2)**2
    end Subroutine jac
    end Module m_newton


Program main
use m_newton 
call solve()
end Program main


微信截图_20231130103719.png (58.81 KB, 下载次数: 121)

这是问题

这是问题
分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

167

帖子

2

主题

1

精华

大师

Vim

F 币
1011 元
贡献
484 点

规矩勋章

沙发
发表于 2023-11-30 16:45:30 | 只看该作者
41行
[Fortran] 纯文本查看 复制代码
Integer::i,j,k,N

2022

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1598 元
贡献
689 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

板凳
发表于 2023-12-1 08:25:13 | 只看该作者
1. 建议每个 module 下面都写 implicit none
2. 建议不要写任何   Implicit real*8(a-z) ,也就是说 implicit 后面只写 none
3. 对contains在module下的过程(函数和子程序),不写 implicit(它会继承module的implicit none)。
4. 没有contains在module和program下面的过程(外部过程),写 implicit none



22

帖子

11

主题

0

精华

入门

F 币
97 元
贡献
54 点
地板
 楼主| 发表于 2023-12-6 15:47:31 | 只看该作者
fcode 发表于 2023-12-1 08:25
1. 建议每个 module 下面都写 implicit none
2. 建议不要写任何   Implicit real*8(a-z) ,也就是说 implic ...

谢谢您的建议

22

帖子

11

主题

0

精华

入门

F 币
97 元
贡献
54 点
5#
 楼主| 发表于 2023-12-6 15:48:04 | 只看该作者

谢谢,已解决
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-11-23 02:09

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表