数值计算编程

第一次上机!Lagrange插值!参数：n插值多项式的阶数，(x,y)插值节点，xx插值点functionLagrange(n,x,y,xx)implicitnoneinteger::n,i,jreal(8),intent(in)::x(0:n),y(0:n),xxreal(8)::Lagrange,lLagrange=0doi=0,nl=y(i)doj=0,nif(j.ne.i)thenl=l*(xx-x(j))/(x(i)-x(j))endifenddoLagrange=Lagrange+lenddoendfunctionLagrange!牛顿插值functionNewton(n,x,y,xx)implicitnoneinteger::n,i,jreal(8),intent(in)::x(0:n),y(0:n),xxreal(8)Newton,P(0:n,0:n),Qp(:,0)=yNewton=p(0,0)doi=1,ndoj=i,np(j,i)=(p(j,i-1)-p(j-1,i-1))/(x(j)-x(j-1))enddoQ=1doj=1,iQ=Q*(xx-x(j-1))enddoNewton=Newton+Q*p(i,i)enddoendfunctionNewton第二次上机realfunctionf(x)implicitnonerealxif(x.eq.0)thenf=1elsef=sin(x)/xendifendfunctionrealfunctiong(x)implicitnonerealxg=4/(1+x**2)endfunctiong!复合梯形求积realfunctionT1(f,a,b,n)implicitnonereal,external::frealx,a,bintegeri,nT1=f(a)+f(b)doi=1,n-1x=a+(b-a)*i/nT1=T1+2*f(x)enddoT1=T1*(b-a)/n/2endfunctionT1!梯形求积逐次对分递推recursiverealfunctionT2(f,a,b,k)result(T2result)implicitnonereal,external::frealy,a,bintegeri,j,kif(k.eq.0)thenT2result=(f(a)+f(b))*(b-a)/2elsey=0doj=1,2**(k-1)y=y+f(a+(2*j-1)*(b-a)/2**k)enddoT2result=T2(f,a,b,k)/2+y*(b-a)/2**kendifendfunctionT2!SimpsonrealfunctionSimpson(f,a,b,k)implicitnonereal,external::freala,b,T2integerkif(k.eq.1)thenSimpson=(4*T2(f,a,b,k)-T2(f,a,b,k-1))/3elseprint*,errorstopendifendfunctionSimpson!CotesrealfunctionCotes(f,a,b,k)realexternal::freala,b,Simpsonintegerkif(k.eq.2)thenCotes=(16*Simpson(f,a,b,k)-Simpson(f,a,b,k-1))/15elseprint*,errorstopendifendfunctionCotes!RombergrealfunctionRomberg(f,a,b,k)realexternal::freala,b,Cotesintegerkif(k.eq.3)thenCotes=(64*Cotes(f,a,b,k)-Cotes(f,a,b,k-1))/63elseprint*,errorstopendifendfunctionRomberg!Romberg求积递推!参数列表:f被积函数，a积分下限，b积分上限，k递推次数，m加速次数recursiverealfunctionR(f,a,b,k,m)result(RR)implicitnonerealexternal::freala,b,yintegerk,m,jif(k.gt.0)thenif(m.gt.0)thenRR=(4**m*R(f,a,b,k,m-1)-R(f,a,b,k-1,m-1))/(4**m-1)elseif(m.eq.0)thenif(k.eq.0)thenRR=(f(a+f(b)))*(b-a)/2elsey=0doj=i,2**(k-1)y=y+f(a+(2*j-1)*(b-a)/2**k)enddoRR=R(f,a,b,k-1,m)/2+(b-a)*y/2**kendifelseprint*,errorstopendifelseprint*,errorstopendifendfunctionR第三次上机!待解方程realfunctionf(x)realxf=x**3-x-1endfunctionf!二分法realfunctionDichotomy(f,limit,a,b)implicitnonereal,external::freala,b,limit,xa,xbif(sign(1.0,f(a))*sign(1.0,f(b)).gt.0)thenprint*,'error:f(a)*f(b)isgreaterthanzero.'endifxa=axb=bdowhile(abs(xa-xb).ge.limit)Dichotomy=(xa+xb)/2if(sign(1.0,f(Dichotomy))*sign(1.0,f(xa)).lt.0)thenxb=Dichotomyelsexa=DichotomyendifenddoDichotomy=(xa+xb)/2endfunctionDichotomy!弦截法realfunctionSecant(f,limit,a,b)implicitnonereal,external::freallimit,a,b,e,x0,x1,x2x0=ax1=be=limitdowhile(e.ge.limit)x2=x1-f(x1)*(x1-x0)/(f(x1)-f(x0))e=abs(x2-x1)x0=x1x1=x2enddoSecant=x2endfunctionSecant