# Demo entry 6645676

f90

Submitted by anonymous on Oct 11, 2017 at 08:58
Language: Fortran. Code size: 4.8 kB.

```program oval1
implicit none
dimension x(10),y(10),B(5),BB(5),XS(5,5),js(5),pfh(10),result(10),o(10),p(10),q(10)
real(kind=8) :: x,y,B,BB,XS,pfh,dx,dy,result,sum_pfh,o,p,q
real(kind=8),parameter :: esp=0.00000001
integer :: i,t,test=0,js

data x/1.02,0.95,0.87,0.77,0.67,0.56,0.44,0.30,0.16,0.01/
data y/0.39,0.32,0.27,0.22,0.18,0.15,0.13,0.12,0.13,0.15/

call jisuanjuzhen(x,y,xs,bb)

call agaus(xs,bb,5,b,test,js)

if(test.ne.0) then
open(unit=9,file="2-1.txt")
call xieru(x,y,t,b,9)
close(9)
endif

write(*,*) "~~~~~~~~~~~~~~~~下面计算误差平方和~~~~~~~~~~~~~~~~~~"
call wucha_pfh(x,y,b,result,pfh,sum_pfh)

stop
end program oval1

subroutine jisuanjuzhen(x,y,xs,bb) !XS点乘B=BB
real(kind=8) x(10),y(10),xs(5,5),bb(5)
real(kind=8) :: x_=0,y_=0,x2=0,x2y=0,x3y=0,xy=0,xy2=0,xy3=0,x2y2=0,y3=0,y2=0,y4=0,x3=0
integer :: i,j,t

do t=1,10,1
x_=x_+x(t)  					!x
y_=y_+y(t) 					!y
x2=x2+x(t)**2 				!x2
y2=y2+y(t)**2 				!y2
x3=x3+x(t)**3 				!x3
y3=y3+y(t)**3 				!y3
xy=xy+x(t)*y(t) 				!xy
x2y=x2y+x(t)*x(t)*y(t)			!x2y
xy2=xy2+x(t)*y(t)*y(t) 			!xy2
x3y=x3y+x(t)**3*y(t)			!x3y
xy3=xy3+x(t)*y(t)**3 		 	!xy3
x2y2=x2y2+x(t)**2*y(t)**2 		!x2y2
y4=y4+y(t)**4 				!y4
enddo

70 format(13f10.6)

XS(1,1)=10   						!给系数矩阵赋值
XS(1,2)=x_
XS(1,3)=y_
XS(1,4)=xy
XS(1,5)=y2

XS(2,2)=x2
XS(2,3)=xy
XS(2,4)=x2y
XS(2,5)=xy2

xs(3,3)=y2
XS(3,4)=xy2
xs(3,5)=y3

xs(4,4)=x2y2
xs(4,5)=xy3

xs(5,5)=y4

do i=2,5
do j=1,i-1
xs(i,j)=xs(j,i)
enddo
enddo

write(*,*) "xs系数矩阵如下："
do i=1,5
write(*,'(5f12.6)') xs(i,:)
enddo
write(*,*)

bb(1)=x2
bb(2)=x3
bb(3)=x2y
bb(4)=x3y
bb(5)=x2y2
end subroutine jisuanjuzhen

subroutine xieru(x,y,t,b,unit)
real(kind=8) x(10),y(10),b(5)
integer :: t,unit

write(unit,*)"x,y:"
do t=1,10
write(unit,"(2f10.6)") x(t),y(t)
enddo 								!在文件中输入X，Y坐标
write(unit,*) 						!换行

do t=1,5
write(unit,20) t-1,B(t)
write(*,20) t-1,B(t)
enddo

write(*,*)
write(*,*)"椭圆轨道公式为："
write(*,10) (b(i),i=1,5)
write(*,*)

write(unit,*)
write(unit,*)"椭圆轨道公式为："
write(unit,10) (b(i),i=1,5)

10 format(f10.6,'+',f10.6,'x+',f10.6,'y+',f10.6,'xy+',f10.6,'y^2=x^2')
20 format(1x,'b',i1,'=',f10.6)
end subroutine xieru

subroutine wucha_pfh(x,y,b,result,pfh,sum_pfh) 							!子程序，求误差平方和
real(kind=8) :: x(10),y(10),result(10),pfh(10),o(10),p(10),q(10),b(5),sum_pfh
integer :: i,j
call qiuxishu(x,b,o,p,q) 					!调用求系数的子程序，o,p,q已经求完
write(*,"('系数',11x,'o',10x,'p',14x,'q')")
do i=1,10
write(*,"(i2,3f14.6)") i,o(i),p(i),q(i) 		!写下o,p,q验证
enddo
write(*,*) "~~~~~~~~~~~~对比y(i),result(i)~~~~~~~~~~~~~~~"
write(*,"('i',10x,'y(i)',7x,'result(i)',6x,'pfh(i)')")
do i=1,10
call solve_eqation(o(i),p(i),q(i),y(i),result(i))	!调用子程序solve_eqation，
pfh(i)=(y(i)-result(i))**2 								!获取误差平方和的值
write(*,"(i2,2f14.6,f14.9)") i,y(i),result(i),pfh(i)
enddo

sum_pfh=sum(pfh)
write(*,*) "~~~~~~~~~~~~~~~~~~误差平方和：~~~~~~~~~~~~~~~~"
write(*,*) sum_pfh

end subroutine wucha_pfh

subroutine qiuxishu(x,b,o,p,q)
real(kind=8) :: x(10),o(10),p(10),q(10),b(5)
integer :: i

do i=1,10
o(i)=b(5)
p(i)=b(3)+b(4)*x(i)
q(i)=b(1)+b(2)*x(i)-x(i)**2
enddo
end subroutine qiuxishu

subroutine solve_eqation(a,b,c,y0,y_nihe) 	!ax2+bx+c=0,y0是题目给的y,y_nihe是返回值
real(kind=8) :: a,b,c,y0,y_nihe,y1,y2,w,u 	!y1,y2是一元二次方程的两个根,w=b2-4ac,u=sqrt(w)

w=b*b-4*a*c
u=sqrt(w)
y1=(-b+u)/(2.0*a)
y2=(-b-u)/(2.0*a)

if(abs(y1-y0) .lt. abs(y2-y0)) then
y_nihe=y1
else
y_nihe=y2
endif
end subroutine solve_eqation

subroutine agaus(a,b,n,x,l,js) 			!求解线代方程
dimension a(n,n),x(n),b(n),js(n)
real(kind=8) :: a,b,x,t
integer :: js
l=1
do 50 k=1,n-1
d=0.0
do 210 i=k,n
do 210 j=k,n
if (abs(a(i,j)).gt.d) then
d=abs(a(i,j))
js(k)=j
is=i
end if
210	  continue
if (d+1.0.eq.1.0) then
l=0
else
if (js(k).ne.k) then
do 220 i=1,n
t=a(i,k)
a(i,k)=a(i,js(k))
a(i,js(k))=t
220	      continue
end if
if (is.ne.k) then
do 230 j=k,n
t=a(k,j)
a(k,j)=a(is,j)
a(is,j)=t
230	      continue
t=b(k)
b(k)=b(is)
b(is)=t
end if
end if
if (l.eq.0) then
write(*,100)
return
end if
do 10 j=k+1,n
a(k,j)=a(k,j)/a(k,k)
10	  continue
b(k)=b(k)/a(k,k)
do 30 i=k+1,n
do 20 j=k+1,n
a(i,j)=a(i,j)-a(i,k)*a(k,j)
20	    continue
b(i)=b(i)-a(i,k)*b(k)
30	  continue
50	continue
if (abs(a(n,n))+1.0.eq.1.0) then
l=0
write(*,100)
return
end if
x(n)=b(n)/a(n,n)
do 70 i=n-1,1,-1
t=0.0
do 60 j=i+1,n
t=t+a(i,j)*x(j)
60	  continue
x(i)=b(i)-t
70	continue
100	format(1x,' fail ')
js(n)=n
do 150 k=n,1,-1
if (js(k).ne.k) then
t=x(k)
x(k)=x(js(k))
x(js(k))=t
end if
150	continue
return

end subroutine agaus
```

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.