最小二乘法-Fortran版本
! https://www.360docs.net/doc/e410551889.html,
subroutine Lsqnonneg(A,b,x,m,n)
!https://www.360docs.net/doc/e410551889.html,
!代码文件调用方式：
!call Lsqnonneg(A,b,x,m,n)
!求解方程组 Ax=b
!其中系数矩阵 A 为 m×n 的矩阵，
!本算法对于 m> n时（即超定方程组）有效。!m !***************** ! 用最小二乘法计算方程组 ! 传进来 Amatrix Bvector use constants use vars_fangchengzu_shuju !use qiujie_jieguo implicit none integer, intent(in) :: m,n logical, allocatable :: p(:),r(:) real, allocatable :: w(:),s(:) real,intent(in),dimension(m,n) :: A real,intent(in),dimension(m) :: b real,intent(inout),dimension(n) :: x real, allocatable :: AT(:,:),bAx(:) ! transpose of A ! b-Ax real, allocatable :: norm_lie(:) ! 放矩阵A中，n列数据的绝对值的和 real :: tol real :: arg_min ! arg_mig=||A^T(b-Ax)|| real :: max_wi integer :: i,jc,j(1),iter integer :: itmax iter=0; !ishape=shape(Amatrix) !m=ishape(1);n=ishape(2) allocate(p(n),r(n)) p=.false. r=.true. allocate(w(n),s(n)) w=0.;s=0. allocate(at(n,m),bAx(m)) !A=Amatrix;b=Bvector; at=transpose(A); !deallocate(Amatrix,Bvector) x=0. !******************** ! 获得误差的判定值 allocate(norm_lie(n)) norm_lie=0. do i=1,n do jc=1,m norm_lie(i)=norm_lie(i)+abs(A(jc,i)); enddo enddo tol=1.0e-9*maxval(norm_lie)*m itmax=3*n; ! 最大的迭代次数。 !****************** bAx=b-matmul(A,x) w=matmul(at,bAx) !**************************************************** do while (any(R) .and. maxval(w)>tol) !***! do while j=MAXLOC(w) ! j=argmax(wi) p(j(1))=.true.; r(j(1))=.false. !****************** iter=iter+1; call lsq_sp(A,b,x,p,r,m,n,tol) ! 返回了新的 x p, r !****************** !****************** bAx=b-matmul(A,x) w=matmul(at,bAx) !******************* write(*,*) "iter=",iter write(*,*) "max_w=",maxval(w) if(iter>itmax) then write(*,*)"迭代次数过多，请增大 tol 的值，即增大误差。" write(*,*)" tol= ",tol return endif enddo !***! do while !**************************************************** return end subroutine Lsqnonneg !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !*****************************************************! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine Lsq_sp(A,b,x,p,r,m,n,tol) ! Lsqnonneg 的子程序 !use qiujie_jieguo !x(n） use imsl implicit none real :: alpha integer :: m,n,i,num_pt real,intent(in) :: tol; integer :: j(1); real,intent(in), dimension(m,n) :: A real,intent(in), dimension(m) :: b real,intent(inout), dimension(n) :: x logical,intent(inout), dimension(n) :: p ,r real, allocatable :: at(:,:) real, allocatable :: s(:) real, allocatable :: sp(:),ap(:,:),apt(:,:),xs(:)!,inv_aa(:,:) ! 其中 xs的维数和s的第二维，也就是p为 true 的个数一致。 ! xs 作为记录 x/(x-s) 的数组，然后取最大值并放入 alpha real,allocatable :: w(:) allocate(w(n),at(n,m)) w=0.; allocate(s(n)) s=0. !************ num_pt=0; do i=1,n if(p(i)==.true.) then num_pt=num_pt+1; endif enddo allocate(sp(num_pt),ap(m,num_pt),apt(num_pt,m), xs(num_pt)) sp=0.;ap=0.;apt=0.;xs=0.;!inv_aa=0. !************ ! 下面计算 sp num_pt=0; do i=1,n if(p(i)==.true.) then num_pt=num_pt+1; ap(:,num_pt)=a(:,i) endif enddo apt=transpose(ap); !inv_aa=.i. (matmul(apt,ap)) sp=matmul(matmul(.i.(matmul(apt,ap)),apt),b) deallocate(apt,ap) !************ num_pt=0; do i=1,n if(p(i)==.true.) then num_pt=num_pt+1; s(i)=sp(num_pt) endif enddo !************ if (minval(sp)<=0.)then num_pt=0; do i=1,n ! write(*,*) "num_pt=",num_pt if(p(i)==.true.) then num_pt=num_pt+1 xs(num_pt)=x(i)/(x(i)-s(i)) endif enddo alpha=-minval(xs); x=x+alpha*(s-x); !****************** !bAx=b-matmul(A,x) w=0. w=matmul(at,(b-matmul(A,x))) if(any(R) .and. maxval(w)>tol) then j=MAXLOC(w) p(j(1))=.true.; r(j(1))=.false. endif ! 更新 RP !****************** endif deallocate(sp,xs) !************ !************ ! 下面再次计算 sp num_pt=0; do i=1,n if(p(i)==.true.) then num_pt=num_pt+1; endif enddo allocate(sp(num_pt),ap(m,num_pt),apt(num_pt,m), xs(num_pt)) sp=0.;ap=0.;apt=0.;xs=0.;!inv_aa=0. !************ num_pt=0; do i=1,n if(p(i)==.true.) then num_pt=num_pt+1; ap(:,num_pt)=a(:,i) endif enddo apt=transpose(ap); !inv_aa=.i. (matmul(apt,ap)) sp=matmul(matmul(.i.(matmul(apt,ap)),apt),b) !************ num_pt=0; s=0. do i=1,n if(p(i)==.true.) then num_pt=num_pt+1; s(i)=sp(num_pt) endif enddo x=s; deallocate(sp,ap,apt, xs,s) return end subroutine Lsq_sp FORTRAN调试程序的时候注意的问题 调试程序的时候注意的问题。 程序编好，能够直接运行而且完全正确，基本不可能，这就有调试程序的问题。主要有一下几个方面： 其一，在每个子程序被调用的时候留个心眼，写个输出语句，表示程序已经运行到了这里。这样一个小提示会给调试带来巨大的方便，如果程序运行出错，至少你可以知道它是在运行到哪里出的错，这样，直接去检查那个程序就可以了。不必重头开始检查。 其二，注意对中间计算结果的输出。有时候，而且是很多的时候，程序编译成功，运行也没有问题，就是结果不对，这肯定是计算原理有问题，此时，输入一些重要步骤的中间结果，往往可以检查出问题所在。同时，就算查出了问题所在，也可以不删除这段输出中间计算结果的代码，有可能后面还会有用处，此时，在每行输出语句前加一个感叹号，把代码变成注释的绿体字就可以了。等到再次需要输出的时候，直接删除一个“！”比再写一遍输出代码，当然要简单的多。 其三，对WATCH功能的应用，FORTRAN提供的这个功能很实用，可以查很多问题，尤其是程序中间计算值，这个和上述的中间结果的输出有点相似。但两者的不同是前者可以进行中间结果的输出控制，就是只有符合了某个条件的才能被输出，这样可以便捷程序的调试，同时对中间结果输出后可以直接用STOP停止程序的运行，这样对于大型程序来说，节省了很多后面继续计算的时间——因为前面的结果已经不对了，后面的算也是白算。 其四，对中间结算结果输出形式的控制，一般来说，FORTRAN计算结果可以输出到文件里面和计算界面两个地方。对较大的计算结果，可以输出到文件里面，反之较少的结果可以直接输出到屏幕上，为了增强数据的可读性，最好进行有格式的数据输出，以利于相同性质的数据的比较。输出到屏幕上的结果直接用WRITE(6,*)就可以（无格式），对于输出到文件里面的数据，可以省些事情，直接用WRITE(X,*)就可以，其中X是一个任意的正整数，最好大于10，也不用事先对这个X设备进行说明，程序会将结果输出到一个FORT.X的文件里面，例如10，就是FORT.10，此时，用NOTEPAD或者ULTRA-EDIT都可以把它打开——FORT.10实质上就是一个.DAT的文件，你可以把它重命名。 3.对数据计算时的误差控制。 以前觉得小数点后的误差不是那回事，没有太在意，可经过实战，终于明白了小小的误差完全可以改变整个计算的结果。因此，如果程序能够输入结果而不正确时，除了寻找算法的问题，不要忽略的误差。一般认为，FORTRAN的REAL变量小数点后8位数字误差定义已经足够，而事实上,这个精度可能在一些情况下不满足，这个时候，需要用更精确的变量类型——REAL(8),同理，当要判断两个数是否相等的时候，一定要慎用相等判断（.EQ.）这个比较运算符，因为任何数据，别看着在现实中它们一定相等，在程序中就不一定了。一旦经过了计算，就不可避免的产生了舍入误差，对于整数和有限几位循环的有理数都问题不大，可一旦是一个无理数或者无限循环的小数，只有在判断了小数点后的每一位都相等的时候，程序才判断为相等成立。这个相等的标准是非常苛刻的，所以一般情况下，可行的方法是将 最小二乘法及其应用 1． 引言 最小二乘法在19世纪初发明后,很快得到欧洲一些国家的天文学家和测地学家的广泛关注。据不完全统计,自1805年至1864年的60年间,有关最小二乘法的研究论文达256篇,一些百科全书包括1837年出版的大不列颠百科全书第7版,亦收入有关方法的介绍。同时,误差的分布是“正态”的,也立刻得到天文学家的关注及大量经验的支持。如贝塞尔( F. W. Bessel, 1784—1846)对几百颗星球作了三组观测,并比较了按照正态规律在给定范围内的理论误差值和实际值,对比表明它们非常接近一致。拉普拉斯在1810年也给出了正态规律的一个新的理论推导并写入其《分析概论》中。正态分布作为一种统计模型,在19世纪极为流行,一些学者甚至把19世纪的数理统计学称为正态分布的统治时代。在其影响下,最小二乘法也脱出测量数据意义之外而发展成为一个包罗极大,应用及其广泛的统计模型。到20世纪正态小样本理论充分发展后,高斯研究成果的影响更加显著。最小二乘法不仅是19世纪最重要的统计方法,而且还可以称为数理统计学之灵魂。相关回归分析、方差分析和线性模型理论等数理统计学的几大分支都以最小二乘法为理论基础。正如美国统计学家斯蒂格勒( S. M. Stigler)所说,“最小二乘法之于数理统计学犹如微积分之于数学”。最小二乘法是参数回归的最基本得方法所以研究最小二乘法原理及其应用对于统计的学习有很重要的意义。 2. 最小二乘法 所谓最小二乘法就是：选择参数10,b b ,使得全部观测的残差平方和最小. 用数学公式表示为： 21022)()(m in i i i i i x b b Y Y Y e --=-=∑∑∑∧ 为了说明这个方法，先解释一下最小二乘原理，以一元线性回归方程为例. i i i x B B Y μ++=10 （一元线性回归方程） 最小二乘法在系统辨识中的应用 王文进 控制科学与控制工程学院 控制理论与控制工程专业 2009010211 摘要：在实际的工程中，经常要对一个系统建立数学模型。很多时候，要面对一个未知的系统，对于这些未知系统，我们所知道的仅仅是它们的一些输入输出数据，我们要根据这些测量的输入输出数据，建立系统的数学模型。由此诞生了系统辨识这门科学，系统辨识就是研究怎样利用对未知系统的输入输出数据建立描述系统的数学模型的科学。系统辨识在工程中的应用非常广泛，系统辨识的方法有很多种，最小二乘法是一种应用及其广泛的系统辨识方法。本文主要讲述了最小二乘估计在系统辨识中的应用。 首先，为了便于介绍，用一个最基本的单输入单输出模型来引入系统辨识中的最小二乘估计。 例如：y = ax + （1） 其中：y、x 可测，为不可测的干扰项，a未知参数。通过N 次实验，得到测量数据y k和x k ，其中k=1、2、3、…，我们所需要做的就是通过这N次实验得到的数据，来确定未知参数a 。在忽略不可测干扰项的前提下，基本的思想就是要使观测点y k和由式（1）确定的估计点y的差的平方和达到最小。用公式表达出来就是要使J最小： 确定未知参数a的具体方法就是令： J a = 0 , 导出 a 通过上面最基本的单输入单输出模型，我们对系统辨识中的最小二乘法有了初步的了解，但在实际的工程中，系统一般为多输入系统，下面就用一个实际的例子来分析。在接下来的表述中，为了便于区分，向量均用带下划线的字母表示。 水泥在凝固过程中，由于发生了一系列的化学反应，会释放出一定的热量。若水泥成分及其组成比例不同，释放的热量也会不同。 水泥凝固放热量与水泥成分的关系模型如下： y = a0+ a1x1+…+ a n x n + 其中，y为水泥凝固时的放热量（卡/克）；x1~x2为水泥的几种成分。 楼主为了减少重复回答问题，特编此帖，并不定期添加和更新内容。 错误难免，欢迎讨论，仅供参考。 很多人问哪里可以找到Fortran编译器，有不少热心学友提供网址，特汇集在这里。虽然俺检验过这些链接，但是它们不一定总有效。 Fortran编译器下载： CVF？ FTN95(License:Freeforpersonaluse) 以下操作，如无特别说明，都是以为例。 1.如何加大Stacksize？ 选Project=>Settings=>Link=>Category:Output=>？ Stackallocations Reserve:这里填新值（默认为1M，若需要10M，则填） 2.如何用Fortran批量生成文件？ 设要生成4000个文件，文件名为AA1-AA4000，如何写循环生成文件,而不用写4000次write 命令呢？ 用内部文件： character(len=80)::filename,form integer::i doi=1,4000 selectcase(i) case(1:9) write(form,'(i1)')i case(10:99) write(form,'(i2)')i case(100:999) write(form,'(i3)')i case(1000:9999) write(form,'(i4)')i endselect write(filename,*)"AA",trim(form),".TXT" open(10,file=filename) write(10,*)i close(10) enddo？ stop end 3.如何用Fortran动态生成输出格式？ 设有一个数组data(100)，输出时，希望每行输出num个数，而num由用户输入，如何实现？ 用内部文件： character(len=80)::form real::data(100) integer::i,num data=(/(i,i=1,100)/)/ read(*,*)num write(form,*)"(",num,"" write(*,form)data stop end 4.MS是不是很垃圾？ 是垃圾，其中Bug太多，多到不可用的地步！ 在这个主题里，换了CVF后问题就没了的人已有相当的数目。 如果你用，遇到莫名其妙的错误，建议换，这是一个比较成熟的编译器。 5.如何用F90/95生成随机数？ 注意： 现在计算机产生的随机数都是伪随机数。 random_number(x)产生一个0到1之间的随机数（x可以是向量），但是每次总是那几个数。用了random_seed()后，系统根据日期和时间随机地提供种子，使得随机数更随机了。 programrandom implicitnone real::x callrandom_seed()!系统根据日期和时间随机地提供种子 callrandom_number(x)!每次的随机数就都不一样了 write(*,*)x stop endprogramrandom 6.函数/子程序超载的例子 一、 最小二乘法（LS ） 辨识系统Z(K+2)=1.5*Z(K+1)-0.7*Z(k)+u(K+1)+0.5*u(k)+v(k) 辨识参数 L T L L T L LS y X X X 1)(-Λ =θ 其中 MAT 程序 >> x=[0 1 0 1 1 0 1 1 1]; >> n=403; >> M=[]; >> for i=1:n temp=xor(x(4),x(9)); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end x(1)=temp; end >> v=randn(1,400); >> z=[]; >> z(1)=-1; >> z(2)=0; >> for i=3:402 z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2); end >> H=zeros(400,4); >> for i=1:400 H(i,1)=-z(i+1); H(i,2)=-z(i); H(i,3)=M(i+1); H(i,4)=M(i); end >> Estimate=inv(H'*H)*H'*(z(3:402))' 辨识参数为： Estimate = -1.4916 1.0364 0.4268 >> 二、最小二乘递推法（RLS） 辨识Z(K+2)=1.5*Z(K+1)-0.7*Z(k)+u(K+1)+0.5*u(k)+v(k) 递推公式： 其中： MATLAB程序： >> x=[0 1 0 1 1 0 1 1 1]; n=403; M=[]; for i=1:n temp=xor(x(4),x(9)); M(i)=x(9); for j=9:-1:2 x(j)=x(j-1); end x(1)=temp; end v=randn(1,400); z=[]; z(1)=-1; z(2)=0; for i=3:402 z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2); end P=100*eye(4); Pstore=zeros(4,401); >> Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)]; >> Theta=zeros(4,401); Theta(:,1)=[3;3;3;3]; >> K=[10;10;10;10]; 移动平均法对小儿科门诊人次的预测分析 目的运用移动平均法对小儿科门诊量进行预测。方法选择我院小儿科2010～2013年各季度门诊量，采用移动平均法建立一元线性回归模型，根据回归模型对2014年各季度数据进行评估预测并做95%置信区间检测。结果2014年各季度小儿科的门诊量均在置信区间内。结论移动平均法可以为我们提供较准确的预测数据，为科室的管理提供决策依据。 标签：移动平均；门诊人次；预测 小儿科门诊是医院里门诊量最大的科室之一，同时也是对季节波动较为敏感的科室，它随着季节变化而呈现有规律性的起伏波动[1]，现根据这种规律性的波动采用移动平均法对门诊人次进行预测，①可以对医院、科室的发展规划提供参考，②可据此来调整门诊医师的做诊班次，在患者高峰来临之前做好准备，减少患者的等待时间，提高患者满意度，因而具有十分重要的意义[2]。 1 资料与方法 1.1一般资料根据我院2010～2013年每季度小儿科门诊人次数据为资料进行统计分析。数据来源我院信息科统计报表，真实可靠。 1.2统计学方法数据采用spss 19.0进行分析，以时间序列的实际值为依据，计算出移动平均的季节比率、校正系数和预测误差，求出长期趋势线，以95%的置信区间对长期趋势值进行预测。 2 结果 2.1计算移动平均比率依据门诊人次（Yt）计算4个季度的移动平均数，再对4项移动平均进行2项移动平均，以剔除时间数列中季节的变动和不规则变动的影响，计算移动平均比率，并计算平均指数的均值，本例为1.0010541；按平均1来校正移动平均比率，本例中校正系数=1/1.0010541。季节修正指数=各季节平均指数×校正系数，并求得调和时间序列（Tt），见表1、表2。 2.2计算调和时间序列并进行季节比率预测调和时间序列等于每一实际门诊数值/季节修正指数。以时间序列号为自变量，以调和时间序列值为应变量建立长期时间趋势序列直线：T=a+b×t，根據最小二乘法原则，采用spss 19.0软件计算：Tt=14984.422+524.141×t。 2.3计算误差及置信区间根据所得直线方程求出各期预测值：Yt=（14984.422+524.141×t）×季节修正指数。并再次利用最小二乘原则计算误差平方e2=（Yt-Yt）2，并对其进行合理性分析得到95%相应置信区间，得到其合理性检验并据此对14年各季度门诊人次进行预测和各季度人次分析，见表3。 ''' 最小二乘法拟合函数曲线f(x) 1、拟合多项式为：y=a0+a1*x+a2*x^2+...+ak*x^k 2、求每个点到曲线的距离之和：Loss=∑(yi-(a0+a1*x+a2*x^2+...+ak*x^k))^2 3、最优化Loss函数，即求Loss对a0,a1,...ak的偏导数为0 3.1、数学解法——求解线性方程组： 整理最优化的偏导数矩阵为：X：含有xi的系数矩阵，A：含有ai的系数矩阵，Y：含有yi的系数矩阵 求解：XA=Y中的A矩阵 3.2、迭代解法——梯度下降法： 计算每个系数矩阵A[k]的梯度，并迭代更新A[k]的梯度 A[k]=A[k]-(learn_rate*gradient) ''' import numpy as np import matplotlib.pyplot as plt plt.rcParams['font.sans-serif']=['SimHei'] plt.rcParams['axes.unicode_minus']=False ''' 高斯列主消元算法 ''' #得到增广矩阵 def get_augmented_matrix(matrix,b): row,col=np.shape(matrix) matrix=np.insert(matrix,col,values=b,axis=1) return matrix #取出增广矩阵的系数矩阵（第一列到倒数第二列） def get_matrix(a_matrix): return a_matrix[:,:a_matrix.shape[1]-1] #选列主元，在第k行后的矩阵里，找出最大值和其对应的行号和列号 def get_pos_j_max(matrix,k): max_v=np.max(matrix[k:,:]) pos=np.argwhere(matrix==max_v) i,_=pos[0] return i,max_v #矩阵的第k行后，行交换 def exchange_row(matrix,r1,r2,k): matrix[[r1,r2],k:]=matrix[[r2,r1],k:] return matrix #消元计算(初等变化) def elimination(matrix,k): row,col=np.shape(matrix) for i in range(k+1,row): m_ik=matrix[i][k]/matrix[k][k] matrix[i]=-m_ik*matrix[k]+matrix[i] return matrix #回代求解 def backToSolve(a_matrix): matrix=a_matrix[:,:a_matrix.shape[1]-1]#得到系数矩阵 b_matrix=a_matrix[:,-1]#得到值矩阵 row,col=np.shape(matrix) x=[None]*col#待求解空间X 1、最小二乘法： 1）（用1 T A A 方法计算逆矩阵） #include double mean,temp=0,variance=0; int i; for(i=0;i 最小二乘参数估计 摘要： 最小二乘的一次性完成辨识算法（也称批处理算法），他的特点是直接利用已经获得的所有（一批）观测数据进行运算处理。这种算法在使用时，占用内存大，离线辨识，观测被辨识对象获得的新数据往往是逐次补充到观测数据集合中去的。在应用一次完成算法时，如果要求在每次新增观测数据后，接着就估计出系统模型的参数，则需要每次新增数据后要重新求解矩阵方程()Z l T l l T l ΦΦΦ-∧=1θ。 最小二乘辩识方法在系统辩识领域中先应用上已相当普及，方法上相当完善，可以有效的用于系统的状态估计，参数估计以及自适应控制及其他方面。 关键词： 最小二乘（Least-squares ），系统辨识（System Identification ） 目录： 1.目的 (1) 2.设备 (1) 3引言 (1) 3.1 课题背景 (1) 4数学模型的结构辨识 (2) 5 程序 (3) 5.1 M 序列子函数 ................................................................................. 错误！未定义书签。 5.2主程序............................................................................................... 错误！未定义书签。 6实验结果： ................................................................................................................................... 3 7参考文献： ................................................................................................. 错误！未定义书签。 1.目的 1.1掌握系统辨识的理论、方法及应用 1.2熟练Matlab 下最小二乘法编程 1.3掌握M 序列产生方法 2.设备 PC 机1台（含Matlab 软件） 3引言 3.1 课题背景 最小二乘理论是有高斯（K.F.Gauss ）在1795年提出：“未知量的最大可能值是这样一个数值，它使各次实际观测值和计算值之间的差值的平方乘以度量其精度的数值以后的和最小。”这就是最小二乘法的最早思想。 最小二乘辨识方法提供一个估算方法，使之能得到一个在最小方差意义上与实验数据最 第一章 FORTRAN程序设计基础第15页 1、2 1．简述程序设计的步骤。 “程序设计”：反映了利用计算机解决问题的全过程，通常要经过以下四个基本步骤：(1)分析问题，确定数学模型或方法；(2)设计算法，画出流程图；(3)选择编程工具，编写程序；(4)调试程序，分析输出结果。 2. 什么是算法？它有何特征？如何描述算法？ 解决问题的方法和步骤称为算法。 算法的五个特征：(1) 有穷性。 (2) 确定性。 (3) 有效性。 (4) 要有数据输入。(5) 要有结果输出。 算法的描述有许多方法，常用的有：自然语言、一般流程图、N－S图等。 第二章顺序结构程序设计 第29页 1、2、3、4、5、6、7、8、9 1．简述符号常量与变量的区别？ 符号常量在程序运行过程中其值不能改变。变量在程序运行过程中其值可以改变。 2. 下列符号中为合法的FORTRAN 90标识符的有哪些? (1) A123B (2) M%10 (3) X_C2 (4) 5YZ (5) X+Y (6) F(X) (7) COS(X) (8) A.2 (9) ‘A’ONE (10) U.S.S.R. (11) min*2 (12) PRINT 3. 下列数据中哪一些是合法的FORTRAN常量？ (1) 9,87 (2) .0 (3) 25.82(4) -356231 (5) 3.57*E2 (6) 3.57E2.1 (7) 3.57E+2(8) 3,57E-2 4. 已知A=2,B=3,C=5(REAL)；且I=2,J=3(INTEGER),求下列表达式的值: (1) A*B+C 表达式的值: 11 (2) A*(B+C) 表达式的值: 16 (3) B/C*A 表达式的值: 1.2 (4) B/(C*A) 表达式的值: 0.3 (5) A/I/J 表达式的值: 0.33 (6) I/J/A 表达式的值: 0 (7) A*B**I/A**J*2 表达式的值: 4.5 (8) C+(B/A)**3/B*2. 表达式的值: 7.25 (9) A**B**I 表达式的值: 512 5. 将下列数学表达式写成相应的FORTRAN表达式： (1) 1E-2 （2）(-B+SQRT(B*B-4*A*C)/(2*A) (3) 1+X+X*X/2+X**3/2/3 (4) COS(ATAN((A**3+B**3)**(1.0/3)/(C*C+1))) (5) EXP(A*X**2+B*X+C) (6) COS(X*Y/SQRT(X*X+Y*Y))**3 6. 用FORTRAN语句完成下列操作： (1) 将变量I的值增加1。I=I+1 (2) I的立方加上J,并将结果保存到I中。 I=I**3+J (3) 将E和F中大者存储到G中。G=Max(E,F) (4) 将两位自然数N的个位与十位互换，得到一个新的数存储到M中(不考虑个位为0的情况) M=MOD(N,10)*10+N/10 第三章选择结构程序设计第43页 1、2、3、5、6、7、9 1．分析下列程序运行结果 (1) LOGICAL P INTEGER I,I1,I2,I3 P=.FALSE. READ*,I I1=MOD(I,10) I2=MOD(I/10,10) I3=I/100 科研探索 知识创新 与109 ——科协论坛?2013年第01期（下）—— 基于移动最小二乘法的列车牵引特性曲线拟合 □ 赵笑龙 徐中伟 朱龙 (同济大学 上海 201804) 摘 要：为了精确地得到列车运行时产生的牵引力，需对列车牵引特性曲线进行拟合。详细阐述移动最小二乘法的基本原理，并以6K 型电力机车为例，分别使用分段最小二乘法和移动最小二乘法进行列车牵引特性曲线的拟合，并将拟合的结果进行对比，通过分析误差的大小及曲线的平滑性，表明该方法的优越性和有效性。关键词：移动最小二乘法牵引特性曲线曲线拟合中图分类号：O241.5 文献标识码：A 文章编号：1007-3973（2013）001-109-03 1引言 曲线拟合在工程设计、计算机图形等方面有着广泛的应用。城市轨道列车牵引力的计算一般分为图解法和分析法，分析法是直接根据物理公式进行计算，求出牵引力，而图解法是通过对牵引特性曲线的拟合，来得到牵引力。由于关于列车牵引力计算的大量公式都是通过经验得来的，因此用公式计算法并不比用图解法得到的结果精确度高。所以我们经常通过图解法来得到列车的牵引力。可以通过测量得到一组关 于速度与牵引力的数据点(),i =1,2…n ，由于不知道牵引特性曲线的原型函数f (x )，而且测量数据点也必然会有误差，因此可以根据离散数据点作出拟合曲线。曲线拟合方法中多使用传统的最小二乘法，即求解一个方程组，使得曲线的误差平方和最小。通过求解该方程组，即可得到期望的拟合曲线。但是传统最小二乘法在拟合的过程中有着较大的计算量，并且在数据量比较大的时候，由于形状比较复杂，可能需要分段拟合，这样会造成拟合曲线的不连续。为了克服上述缺点，本文将使用移动最小二乘法(MLS )来拟合列车牵引曲线。2移动最小二乘法(MLS )的逼近原理 本节主要介绍移动最小二乘法(MLS )的基本原理及其数学表达，以任一计算点处拟合系数的计算为例，来解释移动最小二乘法(MLS )的处理方式及思想。 设全局近似函数为(x )，则它在点x 的局部近似为 (1) 其中：(x )(i =1,2,…,m )是m 次完备单项式基函数，是相应的系数，这些系数是空间坐标x 的函数。基函数(x )常取单项式，在一维问题中： 线性基： ()=[1]，=2(2)二次基： ()=[12]，=3 (3) 为了体现不同节点对计算点的影响，引入了权函数，它是一个紧支撑函数，只和以i 为中心的有效区域内的结点值相关，并与节点到计算点的距离有关。其中，高斯型、紧支撑径向基函数、样条函数和指数型是常用的权函数。 在移动最小二乘法（MLS ）的中，的选取使近似函数在计算点x 领域内对待求函数f (x )的局部逼近误差最小。定义 (4) 令J 取最小值，则 (5) 由此可得 (6) 将式(4)写为矩阵形式，可以得到()=1()()( 7) 其中 (8)()=[(-1)(1)，(-2)(2)，…，(-)()(9) =[(1)，(2)， …， ()](10)()=1()，2()，…，()(11)可见，点x 处的系数有两方面决定：一是该节点值， 二是它的权函数。对不同的节点，随着权函数的变化，系数也跟着不停地变化， 相应的近似函数表示为() =( )()=()1()()(o ) (12) 将移动最小二乘法(MLS )和分段最小二乘法的拟合进行对比，可以发现：分段拟合是将节点域分成了若干段，这就等于对分段的节点分别加权；并且在每个分段内，各节点对拟合系数有着相同的影响，即加权函数取为矩形。而在移动最小二乘法(MLS )中，通过引入紧支撑权函数，对每个计算点附近的节点进行加权，计算点即为中心点，权函数的支撑域即为对中心点有影响的点的范围，随着计算点的改变，权函数的支撑域不断地移动；同时，在每个权函数内，各节点对于计算点的影响与权函数取值有关，因此权函数可以发挥一个软分段的作用。 3移动最小二乘法(MLS )的实施方案 在进行移动最小二乘法(MLS )拟合时,根据具体情况的不同，主要要考虑两方面的因素： (1)基函数的选择：基函数可以分为线性基、二次基和高次基。通常，拟合曲线的平滑性会随着基函数次数的增高越来越好，但是，它会造成就计算量的急剧增大，并且还有可能会出现拟合方程组病态的情况。因此，在拟合过程中，最多取三次基。本文将选择二次基进行实验仿真。(2)权函数的选择：在移动二乘法中，权函数起着至关重要的作用。移动二乘法中的权函数需要满足以下条件：1)在支撑域内，( - ) ＞ 0； FAQ之常见错误 2014-02-02 13:45:35 来源：Fcode研讨团队评论：2点击：4419 本文从编译错误，链接错误，运行时错误，计算结果错误等四个方面介绍了常见的错误及解决思路。适合初学者阅读。 首先应该明确：错误有哪几种？我们当前遇到的是何种错误？ 阐述这些问题前，我们先讨论一下常规的应用程序开发的过程： 1>>编写代码，使用一个或多个源代码文件。 2>>对第一步的每一个源代码文件执行编译操作。得到一个或若干个目标代码。 3>>将目标代码，运行时库（Run-time Library）和其他使用到的函数库链接起来。得到一个可执行文件（EXE 或其他） 4>>编写程序的说明书，必要的（输入）数据文件 5>>将上述得到的结果发布给用户。（发布的方式可以是刻录成光盘，销售，放在网站上供别人下载，或者其他） 6>>用户得到程序后，运行，输入数据，得到计算结果。 对于很多 Fortran 程序员来说，可能用户就是自己，也可能仅仅是自己教研室的同事同学。所以第4，5，6步骤很多时候不明显。而如果使用集成开发环境（IDE）进行开发，第1，2，3步骤又可以一键完成。因此，很多初学者就认为，写程序就是：输入代码，运行，得到结果。这样的理解太狭义。 不管我们面对什么使用者来写代码，程序开发应该是上述的过程。我们的编译器，编译环境，也是为这个过程而设计的。 于是，我们将错误分为四种： 一. 编译错误（发生在第2步） 编译错误，一般是源代码书写格式不正确，不符合语法要求。 二. 链接错误（发生在第3步） 链接错误，一般是源代码结构不完整，运行时库或函数库使用不合理。 三. 运行时错误（发生在第6步） 运行时错误，一般是执行代码时，遇到了事先未料及的错误。比如内存不足了，磁盘空间不够了，输入文件格式不对了，输出文件写入失败了等等。 四. 计算结果不符合预期（程序代码不规范，或不符合你的设想） 计算结果不符合预期，可能性就很多了。语法与你的想法不一致，超出函数库的适用范围，执行流程控制不当等等。 这四种错误，其排查难度依次增大。也就是，编译错误最容易排查和修改，而计算结果不正确，最让人头疼。 为明确解释变量和随机误差各产生的效应是多少，统计上把数据点与它在回归直线上相应位置的差异称为残差，把每个残差的平方和称为残差平方和，它表示随机误差的效应。20()n i i S y y == -∑ 设所示直线议程为y ax b =+，最小二乘法就是示使得残差平方和2 1[()]n i i i M y ax b ==-+∑最 小时a 和b 的值。把M 看作a 和b 的函数，通过求多元函数偏导求最小值时的a 和b 。 1(,)2[()]0n a i i i M M a b y ax b a =?==--+=?∑ 1 (,)2[()]0n b i i i M M a b y a x b b =?==--+=?∑ 即：2111n n n i i i i i i i a x b x x y ===+=∑∑∑ 11 n n i i i i a x nb y ==+=∑∑ 化简得：11122 11()n n n i i i i i i i n n i i i i n x y x y a n x x =====-=-∑∑∑∑∑ b y a x =- 求出其取极值时的a 和b 值取可得直线拟合的方程。其C 语言程序如下： #include ********************************************* fortran中批处理命令的实现函数： 利用systemqq命令（需要调用DFLIB 数据库） ********************************************* 例1： USE DFLIB character*100 CMD LOGICAL(4) res CMD="dir/a-d/b/s "//trim(fPath)//" >"//trim(outPut) res=SYSTEMQQ(CMD) 例2： USE DFLIB LOGICAL(4) result result = SYSTEMQQ('copy e:\dir.txt e:\test\dir.txt') !将e:\dir.txt 复制到e:\test\dir.txt文件中。！****************实例3:复制文件************************* 例3： programmain_pro USE DFLIB implicit none integer,parameter::sta_num=123 character(5),dimension(sta_num)::sta_ID character(500)::filein,fileout character(5000)::cmd logical(4)::judge integer::status,is open(1,file='山东.txt',status='old',action='read',iostat=status) read(1,*) do is=1,sta_num read(1,*) sta_ID(is) filein='Z:\data\降水逐小时数据-戴至修\precip_data\'//sta_ID(is)//'_precip.txt' open(2,file= filein,status='old',action='read',iostat=status) if(status/=0) goto 1000 fileout='Z:\data\降水逐小时数据-戴至修\山东省-降水数据\'//sta_ID(is)//'_precip.txt' cmd='copy '//filein//' '//fileout judge=SYSTEMQQ( cmd) 1000 continue enddo end program imsl7.0中用use linear_operators这句话是会出错的，当时intel论坛上也有人问，后来intel 给出了个X64 imsl的补丁，但是32位的没有。 要使用包含在linear_operators这个库中的函数时，要用use+原函数 例如：上面那个例子把use linear_operators改成use operation_xt即可 另外，imsl7.0引用函数和以前版本不一样的 补充一点： linear_operators这个文件是这样子的（一看就知道原因了）： modulelinear_operators usecond_int usedet_int usediag_int usediagonals_int usefft_int useifft_int useeye_int uselin_eig_self_int uselin_sol_self_int usenorm_int useoperation_i useoperation_ix useoperation_t useoperation_h useoperation_tx useoperation_hx useoperation_x useoperation_xi useoperation_xt useoperation_xh useorth_int userand_int userank_int usesvd_int useunit_int useeig_int usechol_int useisnan_int end module 1. 如何加大Stack size? 选Project => Settings => Link => Category: Output => 移动最小二乘法 2.1 移动最小二乘曲线拟合 将拟合函数表述为如下形式： 1 ()()()()()m T i i i f x p x a x p x a x ===∑, (3) 其中a (x )=(a 1(x ), a 2(x ),…, a m (x ))T 为待定系数，p (x )=(p 1(x ), p 2(x ),…, p m (x ))T 为基函数向量，通常需要选择完备多项式基，例如二维情况 线性基 p (x ) = (1, x , y )T (m =3) 二次基 p (x ) = (1, x , y , x 2, xy , y 2)T (m=6) 为了得到较为精确的局部近似值，需使局部近似值f (x i )和节点值y i 之差平方带权最小，因此残差的离散加权L 2范式为： 2 21 1 ()[()]()[()()]n n T i i i i i i i J w x x f x y w x x p x a x y ===--=--∑∑, (4) 其中n 是求解区域内的节点数，f (x )是拟合函数，w (x －x i )是节点x i 的权函数。 权函数应该是非负的，且随着2 i x x -的增加单调递减，权函数还应该具有紧支性，即 在支持域（x 的影响区域）内不等于0，在支持域之外全为0，一般选用圆形作为权函数的 支持域，半径记为r 。常用的权函数是样条函数，记i s x x '=-，s s r ' =，则三次样条函数形式如下： 23 23 214432 4 41 ()4413 32 0 1. s s s s s s s s s ω?-+≤ ?? ?=-+-<≤??>? ?? (5) 要求出待定系数a (x ),先要使J 取得最小值，先将（4）式写成矩阵形式： J = (Pa (x )－Y )T W (x ) (Pa (x )－Y ) 其中 Y = (y 1, y 2,…, y n )T , W (x ) = diag (w 1(x ), w 2(x ),…, w n (x ))，w i (x ) =()i w x x -. 1121112 22212()()()()()()() ()()m m n n m n p x p x p x p x p x p x P p x p x p x ??????=?? ? ??? . 根据最小二乘原理求得待定系数为： a (x ) = A -1(x )B (x )Y 其中A (x ) = P T W (x ) P , B (x ) = P T W (x )。fortran调试经验
最小二乘法及其应用..
最小二乘法在系统辨识中的应用
fortran常见问题解决
系统辨识
移动平均法对小儿科门诊人次的预测分析
最小二乘算法程序Python版
最小二乘法的编程实现
系统辨识最小二乘参数估计matlab
fortran课后习题答案
基于移动最小二乘法的列车牵引特性曲线拟合
fortran常见错误
最小二乘法C语言程序
fortran中批处理实现
fortran问题
移动最小二乘法