多项式拟合.txt

最小二乘法解多项式拟合问题的fortran编写的源程序

doubleprecision a(100,100),c(100,100),d(100),alpha(100) 
doubleprecision gs(100,100),temp,q,delt(100),x(100),y(100) 
integer gsn 
open(100,file='in.txt',status='old') 
read(100,*) n,m 
read(100,*) (x(i),i=1,n),(y(j),j=1,n) 
close(100) 
do 1 i=1,n 
 do 2 j=1,m+1 
910   a(i,j)=x(i)**(j-1) 
2  continue 
1 continue 
open(200,file='process.txt',status='unknown') 
write(200,*) 'a(m,n):    [','m:',n,'           n:',(m+1),']' 
do 3 i=1,n 
 write(200,35) (a(i,j),j=1,(m+1)) 
3 continue 
do 4 i=1,m+1 
 do 7 j=1,m+1 
920   c(i,j)=0 
  do 8 k=1,n 
   c(i,j)=c(i,j)+x(k)**(i+j-2) 
8   continue 
7  continue 
4 continue 
do 9 i=1,m+1 
 d(i)=0 
 do 11 k=1,n 
929   d(i)=d(i)+y(k)*(x(k)**(i-1)) 
11  continue 
9 continue 
write(200,*) 'c(n,n):    [','n:',(m+1),']' 
do 12 i=1,m+1 
 write(200,35) (c(i,j),j=1,m+1) 
12 continue 
write(200,*) 'd(n):    [','n:',']' 
write(200,35) (d(i),i=1,m+1) 
close(200) 
loop=0 
940 gsn=m+1 
open(300,file='equation.txt',status='unknown') 
write(300,*) 'The added matrix is:' 
do 13 i=1,gsn 
 write(300,35) (c(i,j),j=1,gsn),d(i) 
13 continue 
do 14 i=1,gsn 
 do 15 j=1,gsn 
  gs(i,j)=c(i,j) 
15  continue 
950  gs(i,gsn+1)=d(i) 
14 continue 
do 16 k=1,gsn+1 
 do 17 i=k+1,gsn+1 
  if(dabs(gs(k,k)).lt.dabs(gs(i,j))) then 
   do 18 j=1,gsn+1 
    temp=gs(k,j) 
    gs(k,j)=gs(i,j) 
    gs(i,j)=temp 
18       continue 
  endif 
17     continue 
       temp=gs(k,k) 
 do 19 j=K,gsn+1 
  gs(k,j)=gs(k,j)/temp 
19  continue 
       do 21 i=k,gsn+1 
  temp=gs(i,j) 
  do 22 j=k,gsn+1 
   gs(i,j)=gs(i,j)-temp*gs(k,j) 
22    continue 
21  continue 
 write(300,*) '------------------------------' 
 loop=loop+1 
 write(300,*) 'loop:',loop 
 do 23 i=1,m 
  write(300,35) (gs(i,j),j=1,gsn+1) 
23  continue 
16 continue 
alpha(gsn)=gs(gsn,gsn+1) 
do 24 k=gsn-1,1,-1 
 temp=0.0d0 
 do 25 j=k+1,gsn 
  temp=temp+gs(k,j)*alpha(j) 
25  continue 
24 continue 
write(300,*) '------------------------------' 
write(300,35) (alpha(i),i=1,gsn) 
q=0.0d0 
do 26 i=1,n 
 delt(i)=-y(i) 
 do 27 j=1,m+1 
  delt(i)=delt(i)+(x(i)**(j-1))*alpha(j) 
27     continue 
q=q+delt(i)**2 
26 continue 
write(300,*) '------------------------------' 
write(300,*) 'The final result:' 
write(300,31) (x(i),i=1,n) 
write(300,32) (y(i),i=1,n) 
write(300,33) (delt(i),i=1,n) 
write(300,*) 'q:',q 
close(300) 
write(*,*) (alpha(i),i=1,m+1) 
31 format('x(n):   ',100d30.15) 
32 format('y(n):   ',100d30.15) 
33 format('delt(n):',100d30.15) 
35 format(100d30.15) 
end