Fortran Coder

查看: 8059|回复: 8
打印 上一主题 下一主题

[非线性] 最近用fortran编写最小二乘曲线模拟程序,一直出现问题

[复制链接]

10

帖子

2

主题

0

精华

入门

F 币
68 元
贡献
41 点
跳转到指定楼层
楼主
发表于 2017-8-4 16:28:13 | 只看该作者 |只看大图 回帖奖励 |倒序浏览 |阅读模式
[Fortran] 纯文本查看 复制代码
001program Console1
002implicit none
003integer::i,j,n,status,m,f,k
004real*8 Xmin,Ymin,px,px1,average,yy,dt1,dt2,dt3,z,d1,p,c,d2,g,q,dt,err
005real*8,allocatable::X(:),Y(:),d(:),h(:),a(:)
006real s(20),t(20),b(20)
007open(10,file="E:\Fortran\task\2017.8.3\Console1\Console1\11.dat",status="old")
008k=0
009do while(.true.)
010  read(10,*,iostat=status)
011  if(status/=0)exit
012  k=k+1
013end do
014n=k
015allocate(X(n),Y(n))
016rewind(10)
017do j=1,n
018  read(10,*)X(j),Y(j)
019end do
020Xmin=minval(X)
021Ymin=minval(Y)
022close(10)
023if(Ymin<=0.0)then
024  Y=log10(Y-Ymin+1.0)
025else
026    Y=log10(Y)
027end if
028if(Xmin<=0.0)then
029  X=log10(X-Xmin+1.0)
030else
031  X=log10(X)
032end if
033read (*,*)m
034px=(X(1)+X(n))/2.0
035X=X-px
036f=m-1
037px1=((real(n)+1.0)/sum(x**((2.0)*f)))**((1.0)/(2.0*f))
038x=px1*x
039d=y
040loop1:do
041    call hpir1(x,y,a,n,m,dt1,dt2,dt3)
042    write(9,*)(i,a(i),i=1,m)
043    yy=0
044    average=sum(x)/real(n)
045    outter:do j=1,n
046        h(0)=0
047        inner:do i=1,m
048            h(j)=h(j-1)+a(i)*((x(j)-average)**(i-1))
049        end do inner
050        yy=max(yy,abs(y(j)-h(j)))
051    end do outter
052    if(y(j)>h(j)+yy/real(2.0))then
053            y(j)=h(j)+yy/real(2.0)
054    elseif(y(j)<h(j)-yy/real(2.0))then
055            y(j)=h(j)-yy/real(2.0)
056    end if
057    call RMS(i,y,d,n,err)
058    if(err<1.0e-6.or.j==999)then
059        exit loop1
060    else
061        j=j+1
062    end if
063end do loop1
064open(9,file="result01")
065write(9,400)j,m,n,err,dt1,yy
066400 format(3(i3.3,2x),es12.4,2x,es12.4,2x,f12.4)
067end program
068 
069subroutine RMS(i,y,d,n,err)
070implicit none
071integer,intent(in)::n
072real*8,intent(in)::d(i),y(i)
073real*8 err
074integer::i
075err=0
076do i=1,n
077    err=err+(d(i)-y(i))**2
078end do
079err=sqrt(err/real(n))
080return
081end
082 
083subroutine hpir1(x,y,a,n,m,dt1,dt2,dt3)
084real*8 x(n),y(n),a(m),s(20),t(20),b(20)
085double precision dt1,dt2,dt3,z,d1,p,c,d2,g,q,dt
086do 5 i=1,m
0875 a(i)=0.0
088if (m>n) m=n
089if (m>20) m=20
090z=0.0
091do 10 i=1,n
09210 z=z+x(i)/n
093b(1)=1.0
094d1=n
095p=0.0
096c=0.0
097do 20 i=1,n
098    p=p+(x(i)-z)
099    c=c+y(i)
10020 continue
101c=c/d1
102p=p/d1
103a(1)=c*b(1)
104if (m>1) then
105t(2)=1.0
106t(1)=-p
107d2=0.0
108c=0.0
109g=0.0
110do 30 i=1,n
111    q=x(i)-z-p
112    d2=d2+q*q
113    c=y(i)*q+c
114    g=(x(i)-z)*q*q+g
11530 continue
116c=c/d2
117p=g/d2
118q=d2/d1
119d1=d2
120a(2)=c*t(2)
121a(1)=c*t(1)+a(1)
122endif
123do 100 j=3,m
124    s(j)=t(j-1)
125    s(j-1)=-p*t(j-1)+t(j-2)
126if (j>=4) then
127do 40  k=j-2,2,-1
12840 s(k)=-p*t(k)+t(k-1)-q*b(k)
129end if
130s(1)=-p*T(1)-q*b(1)
131d2=0.0
132c=0.0
133g=0.0
134do 70 i=1,n
135    q=s(j)
136    do 60 k=j-1,1,-1
13760  q=q*(x(i)-z)+s(k)
138    d2=d2+q*q
139    c=y(i)*q+c
140    g=(x(i)-z)*q*q+g
14170 continue
142c=c/d2
143p=g/d2
144q=d2/d1
145d1=d2
146a(j)=c*s(j)
147t(j)=s(j)
148do 80 k=j-1,1,-1
149    a(k)=c*s(k)+a(k)
150    b(k)=t(k)
151    t(k)=s(k)
15280 continue
153100 continue
154dt1=0.0
155dt2=0.0
156dt3=0.0
157do 120 i=1,n
158    q=a(m)
159    do 110 k=m-1,1,-1
160110 q=q*(x(i)-z)+a(k)
161dt=q-y(i)
162if (abs(dt)>dt3) dt3=abs(dt)
163dt1=dt1+dt*dt
164dt2=dt2+abs(dt)
165120 continue
166return
167end


360反馈意见截图16901228253413.png (23.88 KB, 下载次数: 310)

请大神解决

请大神解决
分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

10

帖子

2

主题

0

精华

入门

F 币
68 元
贡献
41 点
沙发
 楼主| 发表于 2017-8-4 16:31:53 | 只看该作者
楼主是fortran新手,参考的徐士良的算法集,可是不知道为什么会报错

10

帖子

2

主题

0

精华

入门

F 币
68 元
贡献
41 点
板凳
 楼主| 发表于 2017-8-4 17:19:28 | 只看该作者
对了,再补充一下,我用的IVF,请大神帮帮忙,程序有点大,谢谢

2033

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1647 元
贡献
709 点

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

地板
发表于 2017-8-5 04:45:48 | 只看该作者
Array a in Main Program has NOT allocated !!!

10

帖子

2

主题

0

精华

入门

F 币
68 元
贡献
41 点
5#
 楼主| 发表于 2017-8-6 17:12:36 | 只看该作者
fcode 发表于 2017-8-5 04:45
Array a in Main Program has NOT allocated !!!

可是程序第五行已经定义了它是可变数组啊

6

帖子

0

主题

0

精华

入门

腐女

F 币
76 元
贡献
49 点
6#
发表于 2017-8-6 17:19:21 | 只看该作者
第五行,只是定义它是“可以”分配的,然而并没有分配。
应该像 15 行,对 X Y 分配那样,同样的对 a 进行分配。

10

帖子

2

主题

0

精华

入门

F 币
68 元
贡献
41 点
7#
 楼主| 发表于 2017-8-7 16:24:39 | 只看该作者
fcode 发表于 2017-8-5 04:45
Array a in Main Program has NOT allocated !!!

可是我添加了allocate a(m)进去了,还是报错

10

帖子

2

主题

0

精华

入门

F 币
68 元
贡献
41 点
8#
 楼主| 发表于 2017-8-7 16:27:45 | 只看该作者
燕雨蔷薇 发表于 2017-8-6 17:19
第五行,只是定义它是“可以”分配的,然而并没有分配。
应该像 15 行,对 X Y 分配那样,同样的对 a 进行 ...

请问添在什么地方呢,我试了好几个地方,包括添加在读取m后面也报错
错误        3        Compilation Aborted (code 1)        E:\Fortran\task\2017.8.3\Console1\Console1\Console1.f90        1       
错误        1         error #5082: Syntax error, found IDENTIFIER 'A' when expecting one of: (        E:\Fortran\task\2017.8.3\Console1\Console1\Console1.f90        55       
错误        2         error #6724: An allocate/deallocate object must have the ALLOCATABLE or POINTER attribute.   [M]        E:\Fortran\task\2017.8.3\Console1\Console1\Console1.f90        55       

10

帖子

2

主题

0

精华

入门

F 币
68 元
贡献
41 点
9#
 楼主| 发表于 2017-8-7 23:57:50 | 只看该作者
我直接把m赋值了,这样就不会出现问题了,可是出现了数组越界问题
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

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

GMT+8, 2025-1-13 10:53

Powered by Discuz! X3.4

© 2013-2025 Comsenz Inc.

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