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.

Delete this entry (admin only).