数值分析大作业三四五六七完整版

合集下载

北航数值分析大作业3

北航数值分析大作业3

数值分析第三次作业1.设计方案对Fredholm积分方程,用迭代法进行求解:()'(())u x A u x=,其中11(())()(,)()A u x g x K x y u y dy-=-⋅⎰对于公式中的积分部分用数值积分方法。

复化梯形积分法,取2601个节点,取迭代次数上限为50次。

实际计算迭代次数为18次,最后算得误差为r= 0.97E-10。

复化Simpson积分法,取迭代次数上限为50次,取2*41+1,即83个节点时能满足精度要求。

实际计算迭代次数为17次,最后的误差为r= 0.97E-10。

Guass积分法选择的Gauss—Legendre法,取迭代次数上限为50次,直接选择8个节点,满足精度要求。

实际计算迭代次数为24次,最后算得误差为r= 0.87E-10。

2.全部源程序module integralimplicit nonecontains!//////////复化梯形subroutine trapezoid(m)implicit noneinteger :: i,j,k,mreal*8 :: x(m+1),u(m+1)real*8 :: sum,sum1,g,r,hreal*8 :: e=1.0e-10h=2./mdo i=1,m+1x(i)=-1.+(i-1)*hend dou=0.02do k=1,50do i=1,m+1sum1=0.g=dexp(x(i)*4.)+(dexp(x(i)+4.)-dexp(-4.-x(i)))/(x(i)+4.)do j=2,msum1=sum1+dexp(x(i)*x(j))*u(j)end dosum=h/2.*(dexp(x(i)*-1.)*u(1)+dexp(x(i)*1.)*u(m+1)+2*sum1)u(i)=g-sumend dor=h/2.*((dexp(x(1)*4)-u(1))**2+(dexp(x(m+1)*4)-u(m+1))**2) do i=2,mr=r+h*(dexp(x(i)*4)-u(i))**2end doif(dabs(r)<=e) exitend dowrite(*,*) kopen(1,file="trapezoid.txt")do i=1,m+1write(1,'(3(f18.12))') x(i),u(i),dexp(x(i)*4.)end dowrite(1,'(4x,a2,e9.2)') "r=",rclose(1)returnend subroutine trapezoid!///////////复化simpsonsubroutine simpson(m)implicit noneinteger :: i,j,k,mreal*8 :: x(2*m+1),u(2*m+1)real*8 :: sum,sum1,sum2,g,r,hreal*8 :: e=1.0e-10h=2./(2.*m)do i=1,2*m+1x(i)=-1.+(i-1)*hend dou=0.02do k=1,50do i=1,2*m+1sum1=0.sum2=0.g=dexp(x(i)*4.)+(dexp(x(i)+4.)-dexp(-4.-x(i)))/(x(i)+4.)do j=1,msum1=sum1+dexp(x(i)*x(2*j))*u(2*j)end dodo j=1,m-1sum2=sum2+dexp(x(i)*x(2*j+1))*u(2*j+1)sum=h/3.*(dexp(x(i)*-1.)*u(1)+dexp(x(i)*1.)*u(2*m+1)+4*sum1+2*sum2) u(i)=g-sumend dor=h/3.*((dexp(x(1)*4)-u(1))**2+(dexp(x(2*m+1)*4)-u(2*m+1))**2)do i=1,mr=r+4.*h/3.*(dexp(x(2*i)*4)-u(2*i))**2end dodo i=1,m-1r=r+2.*h/3.*(dexp(x(2*i+1)*4)-u(2*i+1))**2end doif(dabs(r)<=e) exitend dowrite(*,*) kopen(2,file="simpson.txt")do i=1,2*m+1write(2,'(3(f18.12))') x(i),u(i),dexp(x(i)*4.)end dowrite(2,'(4x,a2,e9.2)') "r=",rclose(2)returnend subroutine simpson!///////////Gauss_Legendre法subroutine Gaussimplicit noneinteger,parameter :: m=8integer :: i,j,kreal*8 :: x(m),u(m),a(m)real*8 :: sum,g,rreal*8 :: e=1.0e-10data x /-0.9602898565,-0.7966664774,-0.5255324099,-0.1834346425,&0.1834346425,0.5255324099,0.7966664774,0.9602898565/data a /0.1012285363,0.2223810345,0.3137066459,0.3626837834,&0.3626837834,0.3137066459,0.2223810345,0.1012285363/u=0.02do k=1,50do i=1,mg=dexp(x(i)*4.)+(dexp(x(i)+4.)-dexp(-4.-x(i)))/(x(i)+4.)do j=1,msum=sum+dexp(x(i)*x(j))*u(j)*a(j)end dou(i)=g-sumend dor=0.do i=1,mr=r+a(i)*(dexp(x(i)*4)-u(i))**2end doif(dabs(r)<=e) exitend dowrite(*,*) kopen(3,file="Gauss.txt")do i=1,mwrite(3,'(3(f18.12))') x(i),u(i),dexp(x(i)*4.)end dowrite(3,'(4x,a2,e9.2)') "r=",rclose(3)returnend subroutine Gaussend module!//////////主程序program mainuse integralimplicit noneinteger :: code1=2600integer :: code2=41call trapezoid(code1)call simpson(code2)call Gaussend program3.各种积分方法的节点和数值解(由于数据太多,在打印时用了较计算时少的有效数字)复化Simpson法4.各方法所得曲线(由于所取节点太多,且精度高,所以图中很难看出各曲线的区别。

数值分析作业(3,4,5,6)

数值分析作业(3,4,5,6)

数值分析作业(3)思考题:1:(a)仅当系数矩阵是病态或奇异的时候,不选主元的Gauss消元法才会失败。

(b) 系数矩阵是对称正定的线性方程组总是良态的;(c) 两个对称矩阵的乘积依然是对称的;(d) 如果一个矩阵的行列式值很小,则它很接近奇异;(e) 两个上三角矩阵的乘积仍然是上三角矩阵;(f) 一个非奇异上三角矩阵的逆仍然是上三角矩阵;(g) 一个奇异矩阵不可能有LU分解;(h) 奇异矩阵的范数一定是零;(i) 范数为零的矩阵一定是零矩阵;(j)一个非奇异的对称阵,如果不是正定的则不能有Cholesky分解。

2: 全主元Gauss消元法与列主元Gauss消元法的基本区别是什么?它们各有什么优点?3:满足下面是我什么条件,可以判定矩阵接近奇异?(a)矩阵的行列式小;(b)矩阵的条件数小;(c)矩阵的范数小;(d)矩阵的条件数大;(e )矩阵的范数大; (f )矩阵的元素小4: Jacobi 迭代法和Gauss-Seidel 迭代法相比(a) 两种的基本差别是什么?(b) 哪种方法更适合并行计算?(c )哪种方法更节省存储空间?(d )Jacobi 迭代法是否运算速度更快?习题:1.对矩阵2112112112A -⎡⎤⎢⎥--⎢⎥=⎢⎥--⎢⎥-⎣⎦,试求A 的Cholesky 分解。

2. 对矩阵12122211111,222221112A A --⎡⎤⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦证明:(1)求解以1A 为系数矩阵的线性方程组,Jacobi 迭代是收敛的,而Gauss-Seidel 迭代发散;(2)求解以2A 为系数矩阵的线性方程组,Jacobi 迭代是发散的,而Gauss-Seidel 迭代收敛。

3.对矩阵11,1a a A a a a a ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦(1) 参数a 取什么值时,矩阵是正定的?(2) a 取什么值,求解以A 为系数矩阵的线性方程组,Jacobi 迭代是收敛的。

北航数值分析大作业三

北航数值分析大作业三

一、题目:关于x, y, t, u, v, w 的下列方程组0.5cos 2.670.5sin 1.070.5cos 3.740.5sin 0.79t u v w x t u v w y t u v w x t u v w y +++-=⎧⎪+++-=⎪⎨+++-=⎪⎪+++-=⎩1、试用数值方法求出f(x, y)在区域 {(,)|00.8,0.5 1.5}D x y x y =≤≤≤≤上的一个近似表达式,0(,)kr s rsr s p x y cx y ==∑要求(,)p x y 一最小的k 值达到以下的精度10202700((,)(,))10i j i j i j f x y p x y σ-===-≤∑∑其中,0.08,0.50.05i j x i y j ==+。

2、计算****(,),(,)i j i j f x y p x y (i = 1, 2, …,8;j = 1, 2,…,5)的值,以观察(,)p x y 逼近(,)f x y 的效果,其中,*i x =0.1i , *j y =0.5+0.2j 。

说明:1、用迭代方法求解非线性方程组时,要求近似解向量()k x 满足()(1)()12||||/||||10k k k x x x --∞∞-≤2、作二元插值时,要使用分片二次代数插值。

3、要由程序自动确定最小的k 值。

4、打印以下内容:●算法的设计方案。

●全部源程序(要求注明主程序和每个子程序的功能)。

●数表:,,i j x y (,)i j f x y (i = 0,1,2,…,10;j = 0,1,2,…,20)。

●选择过程的,k σ值。

●达到精度要求时的,k σ值以及(,)p x y 中的系数rs c (r = 0,1,…,k;s = 0,1,…,k )。

●数表:**,,i j x y ****(,),(,)i j i j f x y p x y (i = 1, 2, ...,8;j = 1, 2, (5)。

数值分析作业(完整版)

数值分析作业(完整版)

的逆阵 A ,用左除命令 A \ E 检验你的结果。
clc clear close all A=[1 1 1 1 1;1 2 3 4 5;1 3 6 10 15;1 4 10 20 35;1 5 15 35 70]; fprintf('对上述矩阵进行列主元素分解:\n') for i=1:1:r-1 [mx,ro]=max(abs(A(i:r,i))); % 寻找a阵第i列的最大值 [A(i,:),A(ro+i-1,:)]=exchange(A(i,:),A(ro+i-1,:)); % 进行行与行交换 for j=i+1:1:r A(j,:)=A(j,:)-A(j,i)/A(i,i)*A(i,:); end A End %--矩阵A的逆阵 A1=inv(A) %--左除验证 E=eye(5); A2=A\E % 5x5单位阵 % A阵的逆矩阵 % 输出每次交换后的A
第一章
1、计算积分 I n
Code: clc clear close all n=9; %--梯形积分法 x=0:0.01:1; y=(x.^n).*exp(x-1); In = trapz(x,y); In2=vpa(In,6) % 6位有效数字 %--高精度积分法 F = @(x1)(x1.^n).*exp(x1-1); s = quad(F,0,1); s1=vpa(s,6)
0
0, 0, 0, 0, 0 。
T
if abs(er(:,i-1))<=e fprintf('在迭代 %d 次之后,满足精度要求,x向量的值如下:\n',i); fprintf('x1=%.5f, x2=%.5f, x3=%.5f, x4=%.5f, x5=%.5f\n',x(1,i),x(2,i),x(3,i),x(4,i),x(5,i)); break end end %--绘图 figure(1) plot(1:1:i,x(1,:),'b',1:1:i,x(2,:),'k',1:1:i,x(3,:),'g',1:1:i,x(4,:), 'r',1:1:i,x(5,:),'c') legend('x1','x2','x3','x4','x5') grid on title('Jacobi迭代法——x值随迭代次数变化曲线') figure(2) plot(1:1:i-1,er(1,:),'b',1:1:i-1,er(2,:),'k',1:1:i-1,er(3,:),'g',1:1: i-1,er(4,:),'r',1:1:i-1,er(5,:),'c') legend('△x1','△x2','△x3','△x4','△x5') grid on title('Jacobi迭代法——△x值随迭代次数变化曲线') %% fprintf('\n-------------Gauss-Seidel迭代法---------------------\n'); U=-(A1-D); L=-(A2-D); DL_1=inv(D-L); M1=DL_1*U; b2=DL_1*b; x1(:,1)=M1*x0+b2; for j=2:1:100 x1(:,j)=M1*x1(:,j-1)+b2; er1(:,j-1)=x1(:,j)-x1(:,j-1); if abs(er1(:,j-1))<=e fprintf('在迭代 %d 次之后,满足精度要求,x向量的值如下:\n',j); fprintf('x1=%.5f, x2=%.5f, x3=%.5f, x4=%.5f, x5=%.5f\n',x1(1,j),x1(2,j),x1(3,j),x1(4,j),x1(5,j)); break end end %--绘图 figure(3) plot(1:1:j,x1(1,:),'b',1:1:j,x1(2,:),'k',1:1:j,x1(3,:),'g',1:1:j,x1(4 ,:),'r',1:1:j,x1(5,:),'c') legend('x1','x2','x3','x4','x5')

北航数值分析全部三次大作业

北航数值分析全部三次大作业

北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。

我们被要求实现各种常用的线性方程组求解算法,例如高斯消元法、LU分解法和迭代法等。

我首先学习了这些算法的原理和实现方法,并借助Python编程语言编写了这些算法的代码。

在实验中,我们使用了不同规模和条件的线性方程组进行测试,并比较了不同算法的性能和精度。

通过这个作业,我深入了解了线性方程组求解的原理和方法,提高了我的编程和数值计算能力。

第二次大作业是关于数值积分的方法。

数值积分是数值分析中的重要内容,它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。

在这个作业中,我们需要实现不同的数值积分算法,例如矩形法、梯形法和辛普森法等。

我学习了这些算法的原理和实现方法,并使用Python编写了它们的代码。

在实验中,我们计算了不同函数的积分值,并对比了不同算法的精度和效率。

通过这个作业,我深入了解了数值积分的原理和方法,提高了我的编程和数学建模能力。

第三次大作业是关于常微分方程的数值解法。

常微分方程是数值分析中的核心内容之一,它可以用于描述众多物理、化学和生物现象。

在这个作业中,我们需要实现不同的常微分方程求解算法,例如欧拉法、龙格-库塔法和Adams法等。

我学习了这些算法的原理和实现方法,并使用Python编写了它们的代码。

在实验中,我们解决了一些具体的常微分方程问题,并比较了不同算法的精度和效率。

通过这个作业,我深入了解了常微分方程的原理和方法,提高了我的编程和问题求解能力。

总的来说,北航数值分析课程的三次大作业非常有挑战性,但也非常有意义。

通过这些作业,我在数值计算和编程方面得到了很大的提升,也更加深入地了解了数值分析的理论和方法。

虽然这些作业需要大量的时间和精力,但我相信这些努力将会对我未来的学习和工作产生积极的影响。

北航数值分析报告大作业第三题(fortran)

北航数值分析报告大作业第三题(fortran)

北航数值分析报告大作业第三题(fortran)“数值分析“计算实习大作业第三题——SY1415215孔维鹏一、计算说明1、将x i=0.08i,y j=0.5+0.05j分别代入方程组(A.3)得到关于t,u,v,w的的方程组,调用离散牛顿迭代子函数求出与x i,y j对应的t i,u j。

2、调用分片二次代数插值子函数在点(t i,u j)处插值得到z(x i,y j)=f(x i,y j),得到数表(x i,y j,f(x i,y j))。

3、对于k=1,2,3,4?,分别调用最小二乘拟合子函数计算系数矩阵c rs 及误差σ,直到满足精度,即求得最小的k值及系数矩阵c rs。

4、将x i?=0.1i,y j?=0.5+0.2j分别代入方程组(A.3)得到关于t?,u?,v?,w?的的方程组,调用离散牛顿迭代子函数求出与x i?,y j?对应的t i?,u j?,调用分片二次代数插值子函数在点(t i?,u j?)处插值得到z?(x i?,y j?)=f(x i?,y j?);调用步骤3中求得的系数矩阵c rs求得p(x i?,y j?),打印数表(x i?,y j?,f(x i?,y j?),p(x i?,y j?))。

二、源程序(FORTRAN)PROGRAM SY1415215DIMENSIONX(11),Y(21),T(6),U(6),Z(6,6),UX(11,21),TY(11,21),FXY(11,21), C(6,6) DIMENSIONX1(8),Y1(5),FXY1(8,5),PXY1(8,5),UX1(8,5),TY1(8,5)REAL(8) X,Y,T,U,Z,FXY,UX,TY,C,E,X1,Y1,FXY1,PXY1,UX1,TY1OPEN (1,FILE='第三题计算结果.TXT')DO I=1,11X(I)=0.08*(I-1)ENDDODO I=1,21Y(I)=0.5+0.05*(I-1)ENDDO!*****求解非线性方程组,得到z=f(t,u)的函数*******DO I=1,11DO J=1,21CALL DISNEWTON_NONLINEAR(X(I),Y(J),UX(I,J),TY(I,J)) ENDDO ENDDO!*************分片二次插值得到z=f(x,y)***********DO I=1,11DO J=1,21CALL INTERPOLATION(UX(I,J),TY(I,J),FXY(I,J))ENDDO ENDDOWRITE (1,"('数表(x,y,f(x,y)):')")WRITE (1,"(3X,'X',7X,'Y',10X,'F(X,Y)')")DO I=1,11DO J=1,21WRITE(1,'(1X,F5.2,2X,F5.3,2X,E20.13)') X(I),Y(J),FXY(I,J) ENDDOWRITE (1,"('')")ENDDO!***********最小二乘拟合得到P(x,y)**************N=11M=21WRITE (1,'(" ","K和σ分别为:")')DO K=1,20CALL LSFITTING(X,Y,FXY,C,N,M,K,K,E) WRITE (1,'(I3,2X,E20.13)') K-1,EIF(ETA).OR.(A(L,K)==TA)) THENTA=A(L,K)TL=LDO J=K,NT(K,J)=A(K,J)A(K,J)=A(TL,J)A(TL,J)=T(K,J)ENDDOTB(K)=B(K)B(K)=B(TL)B(TL)=TB(K)ENDIF ENDDODO I=K+1,NM(I,K)=A(I,K)/A(K,K)A(I,K)=0DO J=K+1,NA(I,J)=A(I,J)-M(I,K)*A(K,J) ENDDOB(I)=B(I)-M(I,K)*B(K)ENDDOENDDO!回代过程X(N)=B(N)/A(N,N)DO K=N-1,1,-1S=0.0DO J=K+1,NS=S+A(K,J)*X(J)ENDDOX(K)=(B(K)-S)/A(K,K)ENDDORETURNEND!***********求向量的无穷数************ SUBROUTINE NORM(X,N,A) DIMENSION X(N)REAL(8) X,AA=ABS(X(1))DO I=2,NIF(ABS(X(I))>ABS(X(I-1))) THENA=ABS(X(I)) ENDIFENDDORETURNEND!**************分片二次代数插值************** SUBROUTINE INTERPOLATION(U,V,W) PARAMETER (N=6,M=6)DIMENSION X(N),Y(M),Z(M,N),LK(3),LR(3)REAL(8) X,Y,Z,H,TREAL(8) U,V,W,LK,LR !U,V分别为插值点处的坐标,W为插值结果INTEGER R!**********************数据赋值********************** DATA Y/0.0,0.2,0.4,0.6,0.8,1.0/DATA X/0.0,0.4,0.8,1.2,1.6,2.0/DATA Z/-0.5,-0.42,-0.18,0.22,0.78,1.5,&&-0.34,-0.5,-0.5,-0.34,-0.02,0.46,&&0.14,-0.26,-0.5,-0.58,-0.5,-0.26,&&0.94,0.3,-0.18,-0.5,-0.66,-0.66,&&2.06,1.18,0.46,-0.1,-0.5,-0.74,&&3.5,2.38,1.42,0.62,-0.02,-0.5/H=0.4T=0.2!******************计算K,R************************* IF(UX(N-1)-H/2) THENK=N-1ELSEDO I=3,N-2IF((U>X(I)-H/2).AND.(UY(M-1)-T/2) THENR=M-1 ELSEDO J=3,M-2IF((V>Y(J)-T/2).AND.(VN) P=N IF(P>20) P=20IF(Q>M) Q=MIF(Q>20) Q=20XX=0YY=0D1=NAPX(1)=0.0DO I=1,NAPX(1)=APX(1)+X(I)ENDDOAPX(1)=APX(1)/D1DO J=1,MV(1,J)=0.0DO I=1,NV(1,J)=V(1,J)+Z(I,J)ENDDOV(1,J)=V(1,J)/D1ENDDOIF(P>1) THEND2=0.0APX(2)=0.0DO I=1,NG=X(I)-APX(1)D2=D2+G*GAPX(2)=APX(2)+(X(I)-XX)*G*G ENDDO APX(2)=APX(2)/D2BX(2)=D2/D1DO J=1,MV(2,J)=0.0DO I=1,NG=X(I)-APX(1)V(2,J)=V(2,J)+Z(I,J)*G ENDDOV(2,J)=V(2,J)/D2ENDDOD1=D2ENDIFDO K=3,PD2=0.0APX(K)=0.0DO J=1,MV(K,J)=0.0ENDDODO I=1,NG1=1.0G2=X(I)-APX(1)DO J=3,KG=(X(I)-APX(J-1))*G2-BX(J-1)*G1 G1=G2 G2=GENDDOD2=D2+G*GAPX(K)=APX(K)+X(I)*G*GDO J=1,M V(K,J)=V(K,J)+Z(I,J)*G ENDDOENDDODO J=1,MV(K,J)=V(K,J)/D2ENDDOAPX(K)=APX(K)/D2BX(K)=D2/D1D1=D2ENDDOD1=MAPY(1)=0.0DO I=1,MAPY(1)=APY(1)+Y(I)ENDDOAPY(1)=APY(1)/D1DO J=1,PU(J,1)=0.0DO I=1,MU(J,1)=U(J,1)+V(J,I) ENDDO U(J,1)=U(J,1)/D1ENDDOIF(Q>1)THEND2=0.0APY(2)=0.0DO I=1,MG=Y(I)-APY(1)D2=D2+G*G APY(2)=APY(2)+(Y(I))*G*G ENDDO APY(2)=APY(2)/D2BY(2)=D2/D1DO J=1,PU(J,2)=0.0DO I=1,MG=Y(I)-APY(1)U(J,2)=U(J,2)+V(J,I)*GENDDOU(J,2)=U(J,2)/D2ENDDOD1=D2ENDIFDO K=3,QD2=0.0APY(K)=0.0DO J=1,PU(J,K)=0.0ENDDODO I=1,MG1=1.0G2=Y(I)-APY(1)DO J=3,KG=(Y(I)-APY(J-1))*G2-BY(J-1)*G1 G1=G2 G2=GENDDOD2=D2+G*GAPY(K)=APY(K)+Y(I)*G*G DO J=1,PU(J,K)=U(J,K)+V(J,I)*G ENDDOENDDODO J=1,PU(J,K)=U(J,K)/D2ENDDOAPY(K)=APY(K)/D2BY(K)=D2/D1D1=D2ENDDOV(1,1)=1.0V(2,1)=-APY(1)V(2,2)=1.0DO I=1,PDO J=1,QA(I,J)=0.0ENDDOENDDODO I=3,QV(I,I)=V(I-1,I-1)V(I,I-1)=-APY(I-1)*V(I-1,I-1)+V(I-1,I-2)IF(I>=4) THENDO K=I-2,2,-1V(I,K)=-APY(I-1)*V(I-1,K)+V(I-1,K-1)-BY(I-1)*V(I-2,K) ENDDO ENDIFV(I,1)=-APY(I-1)*V(I-1,1)-BY(I-1)*V(I-2,1)ENDDO DO I=1,PIF(I==1) THENT(1)=1.0T1(1)=1.0ELSEIF(I==2) THENT(1)=-APX(1)T(2)=1.0T2(1)=T(1)T2(2)=T(2)ELSET(I)=T2(I-1)T(I-1)=-APX(I-1)*T2(I-1)+T2(I-2) IF(I>=4) THENDO K=I-2,2,-1T(K)=-APX(I-1)*T2(K)+T2(K-1)-BX(I-1)*T1(K) ENDDOENDIFT(1)=-APX(I-1)*T2(1)-BX(I-1)*T1(1)T2(I)=T(I)DO K=I-1,1,-1T1(K)=T2(K)T2(K)=T(K)ENDDOENDIFDO J=1,QDO K=I,1,-1DO L=J,1,-1A(K,L)=A(K,L)+U(I,J)*T(K)*V(J,L) ENDDOENDDOENDDOENDDODT1=0.0DO I=1,NX1=X(I)DO J=1,MY1=Y(J)X2=1.0DD=0.0DO K=1,PG=A(K,Q)DO KK=Q-1,1,-1G=G*Y1+A(K,KK)ENDDOG=G*X2DD=DD+GX2=X2*X1ENDDODT=DD-Z(I,J)DT1=DT1+DT*DTENDDOENDDORETURNEND三、计算结果数表(x,y,f(x,y)): X Y UX TY F(X,Y) 0.00 0.500 1.345 0.243 0.17E+000.00 0.550 1.322 0.269 0.66E+000.00 0.600 1.299 0.295 0.35E+000.00 0.650 1.277 0.322 0.94E+000.00 0.700 1.255 0.350 0.30E-020.00 0.750 1.235 0.377 -0.87E-010.00 0.800 1.215 0.406 -0.58E+000.00 0.850 1.196 0.434 -0.72E+000.00 0.900 1.177 0.463 -0.54E+000.00 0.950 1.159 0.492 -0.86E+000.00 1.050 1.125 0.550 -0.74E+00 0.00 1.100 1.109 0.580 -0.06E+00 0.00 1.150 1.093 0.609 -0.00E+00 0.00 1.200 1.0790.639 -0.18E+00 0.00 1.250 1.064 0.669 -0.52E+00 0.00 1.3001.050 0.699 -0.19E+00 0.00 1.350 1.037 0.729 -0.48E+00 0.001.400 1.024 0.759 -0.68E+00 0.00 1.450 1.011 0.790 -0.52E+00 0.00 1.500 1.000 0.820 -0.29E+000.08 0.500 1.415 0.228 0.67E+00 0.08 0.550 1.391 0.253 0.08E+00 0.08 0.600 1.368 0.279 0.02E+00 0.08 0.650 1.346 0.306 0.47E+00 0.08 0.700 1.325 0.333 0.57E+00 0.08 0.750 1.304 0.360 0.48E-01 0.08 0.800 1.284 0.388 -0.73E-01 0.08 0.850 1.265 0.416 -0.16E+00 0.08 0.900 1.246 0.444 -0.29E+00 0.08 0.950 1.229 0.473 -0.36E+00 0.08 1.000 1.211 0.502 -0.08E+00 0.08 1.050 1.194 0.531 -0.29E+00 0.08 1.100 1.178 0.560 -0.78E+00 0.08 1.150 1.163 0.589 -0.93E+00 0.08 1.200 1.148 0.619 -0.44E+00 0.08 1.250 1.133 0.649 -0.92E+00 0.08 1.300 1.119 0.679 -0.71E+000.08 1.400 1.093 0.739 -0.37E+00 0.08 1.450 1.080 0.769-0.83E+00 0.08 1.500 1.068 0.799 -0.92E+000.16 0.500 1.483 0.214 0.31E+00 0.16 0.550 1.460 0.239 0.64E+00 0.16 0.600 1.437 0.264 0.91E+00 0.16 0.650 1.414 0.290 0.06E+00 0.16 0.700 1.393 0.316 0.70E+00 0.16 0.750 1.372 0.343 0.59E+00 0.16 0.800 1.352 0.370 0.12E+00 0.16 0.850 1.333 0.398 0.77E-02 0.16 0.900 1.315 0.426 -0.83E-01 0.16 0.950 1.297 0.454-0.58E+00 0.16 1.000 1.279 0.483 -0.20E+00 0.16 1.050 1.2620.512 -0.11E+00 0.16 1.100 1.246 0.541 -0.74E+00 0.16 1.1501.231 0.570 -0.09E+00 0.16 1.200 1.216 0.600 -0.59E+00 0.16 1.250 1.201 0.629 -0.66E+00 0.16 1.300 1.187 0.659 -0.71E+00 0.16 1.350 1.174 0.689 -0.32E+00 0.16 1.400 1.161 0.718-0.56E+00 0.16 1.450 1.148 0.748 -0.31E+00 0.16 1.500 1.136 0.778 -0.75E+000.24 0.500 1.551 0.201 0.66E+01 0.24 0.550 1.527 0.2250.03E+000.24 0.650 1.482 0.275 0.64E+00 0.24 0.700 1.460 0.3010.47E+00 0.24 0.750 1.439 0.327 0.34E+00 0.24 0.800 1.419 0.354 0.24E+00 0.24 0.850 1.400 0.381 0.69E+00 0.24 0.900 1.381 0.409 0.04E-01 0.24 0.950 1.363 0.437 -0.42E-01 0.24 1.000 1.346 0.465 -0.06E+00 0.24 1.050 1.329 0.494 -0.59E+00 0.24 1.100 1.313 0.523 -0.83E+00 0.24 1.150 1.297 0.552 -0.15E+00 0.24 1.200 1.282 0.581 -0.19E+00 0.24 1.250 1.267 0.610 -0.84E+00 0.24 1.300 1.253 0.640 -0.66E+00 0.24 1.350 1.240 0.669 -0.30E+00 0.24 1.400 1.227 0.699 -0.86E+00 0.24 1.450 1.214 0.729 -0.84E+00 0.24 1.500 1.202 0.759 -0.77E+000.32 0.500 1.617 0.188 0.28E+01 0.32 0.550 1.593 0.212 0.49E+01 0.32 0.600 1.570 0.236 0.68E+00 0.32 0.650 1.547 0.261 0.75E+00 0.32 0.700 1.526 0.286 0.60E+00 0.32 0.750 1.505 0.312 0.77E+00 0.32 0.800 1.485 0.339 0.05E+00 0.32 0.850 1.466 0.365 0.99E+00 0.32 0.900 1.447 0.393 0.27E+00 0.32 1.000 1.411 0.448 -0.01E-02 0.32 1.050 1.395 0.477-0.41E-01 0.32 1.100 1.378 0.505 -0.18E+00 0.32 1.150 1.3630.534 -0.25E+00 0.32 1.200 1.347 0.563 -0.29E+00 0.32 1.2501.333 0.592 -0.90E+00 0.32 1.300 1.319 0.621 -0.00E+00 0.32 1.350 1.305 0.650 -0.40E+00 0.32 1.400 1.292 0.680 -0.54E+00 0.32 1.450 1.279 0.710 -0.79E+00 0.32 1.500 1.267 0.739-0.91E+000.40 0.500 1.681 0.177 0.91E+01 0.40 0.550 1.658 0.1990.00E+01 0.40 0.600 1.634 0.223 0.83E+01 0.40 0.650 1.612 0.247 0.02E+01 0.40 0.700 1.591 0.272 0.94E+00 0.40 0.750 1.570 0.298 0.49E+00 0.40 0.800 1.550 0.324 0.94E+00 0.40 0.850 1.530 0.350 0.40E+00 0.40 0.900 1.512 0.377 0.33E+00 0.40 0.950 1.493 0.405 0.99E+00 0.40 1.000 1.476 0.432 0.68E+00 0.40 1.050 1.459 0.460 0.08E-01 0.40 1.100 1.443 0.488 -0.84E-01 0.40 1.150 1.427 0.517-0.98E+00 0.40 1.200 1.412 0.545 -0.27E+00 0.40 1.250 1.397 0.574 -0.06E+000.40 1.350 1.369 0.632 -0.66E+00 0.40 1.400 1.356 0.662-0.37E+00 0.40 1.450 1.343 0.691 -0.43E+00 0.40 1.500 1.331 0.721 -0.12E+000.48 0.500 1.745 0.166 0.69E+01 0.48 0.550 1.721 0.188 0.02E+01 0.48 0.600 1.698 0.211 0.74E+01 0.48 0.650 1.676 0.235 0.40E+01 0.48 0.700 1.654 0.259 0.23E+01 0.48 0.750 1.634 0.284 0.56E+00 0.48 0.800 1.613 0.310 0.28E+00 0.48 0.850 1.594 0.336 0.49E+00 0.48 0.900 1.575 0.363 0.31E+00 0.48 0.950 1.557 0.390 0.66E+00 0.48 1.000 1.539 0.417 0.30E+00 0.48 1.050 1.522 0.444 0.34E+00 0.48 1.100 1.506 0.472 0.07E-01 0.48 1.150 1.490 0.500 -0.62E-01 0.48 1.200 1.475 0.529 -0.45E+00 0.48 1.250 1.460 0.557 -0.86E+00 0.48 1.300 1.446 0.586 -0.39E+00 0.48 1.350 1.432 0.615 -0.22E+00 0.48 1.400 1.419 0.644 -0.67E+00 0.48 1.450 1.406 0.674-0.55E+00 0.48 1.500 1.394 0.703 -0.14E+000.56 0.500 1.808 0.156 0.48E+010.56 0.600 1.761 0.200 0.10E+01 0.56 0.650 1.739 0.2230.68E+01 0.56 0.700 1.717 0.247 0.94E+01 0.56 0.750 1.696 0.272 0.33E+01 0.56 0.800 1.676 0.297 0.11E+00 0.56 0.850 1.657 0.323 0.63E+00 0.56 0.900 1.638 0.349 0.97E+00 0.56 0.950 1.620 0.375 0.52E+00 0.56 1.000 1.602 0.402 0.56E+00 0.56 1.050 1.585 0.429 0.47E+00 0.56 1.100 1.568 0.457 0.20E+00 0.56 1.150 1.552 0.485 0.13E+00 0.56 1.200 1.537 0.513 0.09E-01 0.56 1.250 1.522 0.541 -0.47E-01 0.56 1.300 1.508 0.570 -0.99E+00 0.56 1.350 1.4940.599 -0.82E+00 0.56 1.400 1.481 0.627 -0.26E+00 0.56 1.4501.468 0.657 -0.71E+00 0.56 1.500 1.455 0.686 -0.98E+000.64 0.500 1.870 0.147 0.74E+01 0.64 0.550 1.846 0.1680.10E+01 0.64 0.600 1.823 0.190 0.54E+01 0.64 0.650 1.801 0.213 0.42E+01 0.64 0.700 1.779 0.236 0.56E+01 0.64 0.750 1.758 0.260 0.03E+01 0.64 0.800 1.738 0.285 0.42E+01 0.64 0.850 1.718 0.310 0.41E+010.64 0.950 1.681 0.362 0.36E+00 0.64 1.000 1.664 0.388 0.18E+00 0.64 1.050 1.646 0.415 0.28E+00 0.64 1.100 1.630 0.443 0.07E+00 0.64 1.150 1.614 0.470 0.66E+00 0.64 1.200 1.598 0.498 0.09E+00 0.64 1.250 1.584 0.526 0.50E-01 0.64 1.300 1.569 0.554 -0.88E-01 0.64 1.350 1.555 0.583 -0.76E+00 0.64 1.400 1.542 0.611 -0.66E+00 0.64 1.450 1.529 0.640 -0.33E+00 0.64 1.500 1.516 0.669 -0.56E+00 0.72 0.500 1.931 0.139 0.94E+01 0.72 0.550 1.907 0.159 0.84E+01 0.72 0.600 1.884 0.181 0.36E+01 0.72 0.650 1.862 0.203 0.40E+01 0.72 0.700 1.840 0.226 0.47E+01 0.72 0.750 1.819 0.249 0.56E+01 0.72 0.800 1.799 0.273 0.19E+01 0.72 0.850 1.779 0.298 0.37E+01 0.72 0.900 1.760 0.323 0.86E+01 0.72 0.950 1.742 0.349 0.76E+00 0.72 1.000 1.724 0.375 0.24E+00 0.72 1.050 1.707 0.402 0.55E+00 0.72 1.100 1.691 0.429 0.97E+00 0.72 1.150 1.675 0.456 0.27E+00 0.72 1.200 1.659 0.484 0.31E+000.72 1.300 1.630 0.539 0.49E+00 0.72 1.350 1.616 0.5680.72E-02 0.72 1.400 1.602 0.596 -0.69E-01 0.72 1.450 1.589 0.625 -0.67E+00 0.72 1.500 1.576 0.653 -0.20E+000.80 0.500 1.992 0.131 0.31E+01 0.80 0.550 1.968 0.1510.44E+01 0.80 0.600 1.945 0.172 0.41E+01 0.80 0.650 1.922 0.193 0.45E+01 0.80 0.700 1.900 0.216 0.00E+01 0.80 0.750 1.879 0.239 0.10E+01 0.80 0.800 1.859 0.263 0.16E+01 0.80 0.850 1.840 0.287 0.52E+01 0.80 0.900 1.821 0.312 0.02E+01 0.80 0.950 1.802 0.337 0.38E+01 0.80 1.000 1.784 0.363 0.89E+01 0.80 1.050 1.767 0.389 0.28E+00 0.80 1.100 1.751 0.416 0.09E+00 0.80 1.150 1.734 0.4430.23E+00 0.80 1.200 1.719 0.470 0.93E+00 0.80 1.250 1.704 0.498 0.15E+00 0.80 1.300 1.689 0.525 0.86E+00 0.80 1.350 1.675 0.553 0.64E+00 0.80 1.400 1.662 0.582 0.74E-01 0.80 1.450 1.649 0.610 -0.37E-01 0.80 1.500 1.636 0.638 -0.81E+00K和σ分别为:0 0.93E+031 0.61E+012 0.92E-023 0.53E-034 0.16E-055 0.77E-07系数矩阵Crs(按行)为:0.00E+01 -0.83E+01 0.56E+00 0.97E+00 -0.03E+00 0.70E-010.91E+01 -0.99E+00 -0.96E+01 0.17E+01 -0.66E+00 0.10E-01 0.77E+00 0.42E+01 -0.10E+00 -0.81E+00 0.81E+00 -0.62E-01-0.25E+00 -0.21E+00 0.97E+00 -0.18E+00 0.49E+00 -0.63E-010.34E+00 -0.56E+00 0.69E-01 0.51E+00 -0.77E-01 0.27E-01-0.94E-01 0.94E+00 -0.58E+00 0.69E-01 -0.50E-01 0.53E-02 数表(x,y,f(x,y),p(x,y)):X Y F(X,Y) P(X,Y)0.100 0.700 0.58E+00 0.05E+000.100 1.100 -0.66E+00 -0.26E+00 0.100 1.300 -0.68E+00-0.31E+00 0.100 1.500 -0.52E+00 -0.49E+000.200 0.700 0.54E+00 0.19E+00 0.200 0.900 -0.63E-01 -0.65E-01 0.200 1.100 -0.90E+00 -0.90E+00 0.200 1.300 -0.84E+00 -0.90E+00 0.200 1.500 -0.03E+00 -0.04E+000.300 0.700 0.82E+00 0.09E+00 0.300 0.900 0.48E+00 0.11E+00 0.300 1.100 -0.63E+00 -0.88E+00 0.300 1.300 -0.72E+00 -0.96E+00 0.300 1.500 -0.34E+00 -0.84E+000.400 0.700 0.79E+00 0.89E+00 0.400 0.900 0.56E+00 0.63E+00 0.400 1.100 -0.83E-01 -0.04E-01 0.400 1.300 -0.72E+00 -0.71E+00 0.400 1.500 -0.85E+00 -0.07E+000.500 0.700 0.56E+01 0.92E+01 0.500 0.900 0.51E+00 0.23E+00 0.500 1.100 0.59E+00 0.27E+00 0.500 1.300 -0.53E+00 -0.11E+00 0.500 1.500 -0.67E+00 -0.33E+000.600 0.900 0.14E+00 0.75E+00 0.600 1.100 0.19E+00 0.32E+00 0.600 1.300 -0.70E-01 -0.82E-01 0.600 1.500 -0.08E+00 -0.75E+00 0.700 0.700 0.89E+01 0.29E+01 0.700 0.900 0.91E+01 0.11E+010.700 1.100 0.60E+00 0.97E+00 0.700 1.300 0.22E-01 0.06E-01 0.7001.500 -0.53E+00 -0.80E+00 0.800 0.700 0.09E+01 0.06E+01 0.800 0.900 0.32E+01 0.50E+01 0.800 1.100 0.03E+00 0.79E+00 0.800 1.300 0.25E+00 0.50E+00 0.800 1.500 -0.14E+00 -0.28E+00。

BUAA数值分析大作业三

BUAA数值分析大作业三

北京航空航天大学2020届研究生《数值分析》实验作业第九题院系:xx学院学号:姓名:2020年11月Q9:方程组A.4一、 算法设计方案(一)总体思路1.题目要求∑∑===k i kj s r rsy x cy x p 00),(对f(x, y) 进行拟合,可选用乘积型最小二乘拟合。

),(i i y x 与),(i i y x f 的数表由方程组与表A-1得到。

2.),(**j i y x f 与1使用相同方法求得,),(**j i y x p 由计算得出的p(x,y)直接带入),(**j i y x 求得。

1. ),(i i y x 与),(i i y x f 的数表的获得对区域D ={ (x,y)|1≤x ≤1.24,1.0≤y ≤1.16}上的f (x , y )值可通过xi=1+0.008i ,yj=1+0.008j ,得到),(i i y x 共31×21组。

将每组带入A4方程组,即可获得五个二元函数组,通过简单牛顿迭代法求解这五个二元数组可获得z1~z5有关x,y 的表达式。

再将),(i i y x 分别带入z1~z5表达式即可获得f(x,y)值。

2.乘积型最小二乘曲面拟合2.1使用乘积型最小二乘拟合,根据k 值不用,有基函数矩阵如下:⎪⎪⎪⎭⎫ ⎝⎛=k i i k x x x x B 0000 , ⎪⎪⎪⎭⎫ ⎝⎛=k j jk y y y y G 0000数表矩阵如下:⎪⎪⎪⎭⎫⎝⎛=),(),(),(),(0000j i i j y x f y x f y x f y x f U记C=[rs c ],则系数rs c 的表达式矩阵为:11-)(-=G G UG B B B C T TT )(通过求解如下线性方程,即可得到系数矩阵C 。

UG B G G C B B T T T =)()(2.2计算),(),,(****j i j i y x p y x f (i =1,2,…,31 ; j =1,2,…,21) 的值),(**j i y x f 的计算与),(j i y x f 相同。

数值分析大作业

数值分析大作业

数值分析大作业数值分析大作业姓名:黄晨晨学号:S1*******学院:储运与建筑工程学院学院班级:储建研17-2实验3.1 Gauss消去法的数值稳定性实验实验目的:理解高斯消元过程中出现小主元即很小时引起方程组解数值不定性实验内容:求解方程组Ax=b,其中(1)A1=0.3×10?1559.14315.291?6.130?1211.29521211,b1=59.1746.7812;(2)A2=10?7013 2.099999999999625?15?10102,b2=85.90000000000151;实验要求:(1)计算矩阵的条件数,判断系数矩阵是良态的还是病态的(2)用Gauss列主元消去法求得L和U及解向量x1,x2∈R4(3)用不选主元的高斯消去法求得L和U及解向量x1,x2∈R4(4)观察小主元并分析对计算结果的影响(1)计算矩阵的条件数,判断系数矩阵是良态的还是病态的代码:format longeformat compactA1=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1] b1=[59.17;46.78;1;2]n=4C1=cond(A1,1) %C1为A1矩阵1范数下的条件数C2=cond(A1,2) %C2为A1矩阵2范数下的条件数C3=cond(A1,inf) %C3为1矩阵谱范数下的条件数结果:C1 =1.362944708720448e+02C2 =6.842955771253409e+01C3 =8.431146*********e+01显然A1矩阵为病态矩阵将矩阵A2,b2输入上述代码中求得A2矩阵的条件数为:C1 =1.928316831682894e+01C2 =8.993938090170119e+00C3 =1.835643564356072e+01显然A2矩阵也为病态矩阵(2)用Gauss列主元消去法求得L和U及解向量x1,x2∈R4代码:for k=1:n-1a=max(abs(A1(k:n,k)))if a==0returnendfor i=k:nif abs(A1(i,k))==ay=A1(i,:)A1(i,:)=A1(k,:)A1(k,:)=yx=b1(i,:)b1(i,:)=b1(k,:)b1(k,:)=xbreakendendif A1(k,k)~=0A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k)A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n) elsebreakendendL=tril(A1,0);for i=1:nL(i,i)=1;endLU=triu(A1,0)y1=L\b1x1=U\y1得到如下结果:x1 =3.845714853511634e+001.609517394778522e+00-1.547605454206655e+011.041130489899787e+01将A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2]b2=[8;5.900000000001;5;1]代入上述代码求得结果如下:X2 =4.440892098500626e-16-9.999999999999993e-019.999999999999997e-011.000000000000000e+00(3)用不选主元的高斯消去法求得L和U及解向量x1,x2∈R4代码:format longeformat compactA1=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1] b1=[59.17;46.78;1;2][L,U]=lu(A1)y1=L\b1x1=U\y1求得如下结果:x1=3.845714853511634e+001.609517394778522e+00-1.547605454206655e+011.041130489899787e+01将A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2] b2=[8;5.900000000001;5;1]代入上述代码,求得结果如下:x 2 =4.440892098500626e-16 -9.999999999999993e-01 9.999999999999997e-01 9.999999999999999e-01(2)(3)求得结果相同,可验证结果正确。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

数值分析大作业三四五六七HEN system office room 【HEN16H-HENS2AHENS8Q8-HENH1688】大作业 三1.给定初值0x 及容许误差,编制牛顿法解方程f (x )=0的通用程序.解:Matlab 程序如下:函数m 文件:function Fu=fu(x) Fu=x^3/3-x; end函数m 文件:function Fu=dfu(x) Fu=x^2-1; end用Newton 法求根的通用程序 clear;x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1;while flag==1x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)<ep flag=0; end x0=x1; endfprintf('方程的一个近似解为:%f\n',x0); 寻找最大δ值的程序: cleareps=input('请输入搜索精度:'); ep=input('请输入容许误差:'); flag=1; k=0; x0=0; while flag==1 sigma=k*eps; x0=sigma; k=k+1; m=0; flag1=1;while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0);if abs(x1-x0)<epend m=m+1; x0=x1; endif flag1==1||abs(x0)>=ep flag=0; end endfprintf('最大的sigma 值为:%f\n',sigma);2.求下列方程的非零根5130.6651()ln 05130.665114000.0918x x f x x +⎛⎫=-= ⎪-⨯⎝⎭解:Matlab 程序为:(1)主程序 clear clc format long x0=765; N=100;errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1; while n<Nx=x0-f(x0)/subs(df(),x0); if abs(x-x0)>errorlim n=n+1; else break ; end x0=x; enddisp(['迭代次数: n=',num2str(n)])disp(['所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)]) (2)子函数 非线性函数f function y=f(x)y=log((513+*x)/*x))-x/(1400*; end(3)子函数 非线性函数的一阶导数df function y=df() syms x1y=log((513+*x1)/*x1))-x1/(1400*; y=diff(y);运行结果如下:迭代次数: n=5所求非零根: 正根x1= 负根x2=大作业 四试编写MATLAB 函数实现Newton 插值,要求能输出插值多项式. 对函数21()14f x x =+在区间[-5,5]上实现10次多项式插值.分析:(1)输出插值多项式。

(2)在区间[-5,5]内均匀插入99个节点,计算这些节点上函数f (x )的近似值,并在同一张图上画出原函数和插值多项式的图形。

(3)观察龙格现象,计算插值函数在各节点处的误差,并画出误差图。

解:Matlab 程序代码如下:%此函数实现y=1/(1+4*x^2)的n 次Newton 插值,n 由调用函数时指定 %函数输出为插值结果的系数向量(行向量)和插值多项式 function [t y]=func5(n) x0=linspace(-5,5,n+1)'; y0=1./(1.+4.*x0.^2); b=zeros(1,n+1); for i=1:n+1 s=0; for j=1:i t=1; for k=1:i if k~=jt=(x0(j)-x0(k))*t; end ; end ; s=s+y0(j)/t; end ; b(i)=s; end ;t=linspace(0,0,n+1); for i=1:ns=linspace(0,0,n+1);s(n+1-i:n+1)=b(i+1).*poly(x0(1:i)); t=t+s; end ;t(n+1)=t(n+1)+b(1); y=poly2sym(t);10次插值运行结果: [b Y]=func5(10)Columns 1 through 4Columns 5 through 8Columns 9 through 11Y =b为插值多项式系数向量,Y为插值多项式。

插值近似值:x1=linspace(-5,5,101);x=x1(2:100);y=polyval(b,x)y =Columns 1 through 12Columns 13 through 24Columns 25 through 36Columns 37 through 48Columns 49 through 60Columns 61 through 72Columns 73 through 84Columns 85 through 96Columns 97 through 99绘制原函数和拟合多项式的图形代码:plot(x,1./(1+4.*x.^2))hold allplot(x,y,'r')xlabel('X')ylabel('Y')title('Runge现象')gtext('原函数')gtext('十次牛顿插值多项式')绘制结果:误差计数并绘制误差图:hold offey=1./(1+4.*x.^2)-yey =Columns 13 through 24Columns 25 through 36Columns 37 through 48Columns 49 through 60Columns 61 through 72Columns 73 through 84Columns 85 through 96Columns 97 through 99plot(x,ey)xlabel('X')ylabel('ey')title('Runge现象误差图')输出结果为:大作业五解:Matlab程序为:x = [-520,-280,,-78,,,0,,,78,,280,520]';y = [0,-30,-36,-35,,,0,,,35,36,30,0]';n = 13;%求解Mfor i = 1:1:n-1h(i) = x(i+1)-x(i);endfor i = 2:1:n-1a(i) = h(i-1)/(h(i-1)+h(i));b(i) = 1-a(i);c(i) = 6*((y(i+1)-y(i))/h(i)-(y(i)-y(i-1))/h(i-1))/(h(i-1)+h(i)); enda(n) = h(n-1)/(h(1)+h(n-1));b(n) = h(1)/(h(1)+h(n-1));c(n) = 6/(h(1)+h(n-1))*((y(2)-y(1))/h(1)-(y(n)-y(n-1))/h(n-1));A(1,1) = 2;A(1,2) = b(2);A(1,n-1) = a(2);A(n-1,n-2) = a(n);A(n-1,n-1) = 2;for i = 2:1:n-2A(i,i) = 2;A(i,i+1) = b(i+1);A(i,i-1) = a(i+1);endC = c(2:n);C = C';m = A\C;M(1) = m(n-1);M(2:n) = m;xx = -520::520;for i = 1:51for j = 1:1:n-1if x(j) <= xx(i) && xx(i) < x(j+1)break;endendyy(i) = M(j+1)*(xx(i)-x(j))^3/(6*h(j))-M(j)*(xx(i)-x(j+1))^3/(6*h(j))+(y(j+1)-M(j+1)*h(j)^2/6)*(xx(i)-x(j))/h(j)-(y(j)-M(j)*h(j)^2/6)*(xx(i)-x(j+1))/h(j); end;for i = 52:101yy(i) = -yy(102-i);end;for i = 1:50xx(i) = -xx(i);endplot(xx,yy);hold on;for i = 1:1:n/2x(i) = -x(i);endplot(x,y,'bd');title('机翼外形曲线');输出结果:运行文件,得到2. (1)编制求第一型3 次样条插值函数的通用程序;(2)已知汽车门曲线型值点的数据如下:解:(1)Matlab编制求第一型3 次样条插值函数的通用程序:function [Sx] = Threch(X,Y,dy0,dyn)%X为输入变量x的数值%Y为函数值y的数值%dy0为左端一阶导数值%dyn为右端一阶导数值%Sx为输出的函数表达式n = length(X)-1;h = zeros(1,n-1);f1 = zeros(1,n-1);f2 = zeros(1,n-2);for i = 1:n %求函数的一阶差商h(i) = X(i+1)-X(i);f1(i) = (Y(i+1)-Y(i))/h(i);endfor i = 2:n %求函数的二阶差商f2(i) = (f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i) = 6*f2(i);endd(1) = 6*(f1(1)-dy0)/h(1);d(n+1) = 6*(dyn-f1(n-1))/h(n-1);%赋初值A = zeros(n+1,n+1);B = zeros(1,n-1);C = zeros(1,n-1);for i = 1:n-1B(i) = h(i)/(h(i)+h(i+1));C(i) = 1-B(i);endA(1,2) = 1;A(n+1,n) = 1;for i = 1:n+1A(i,i) = 2;endfor i = 2:nA(i,i-1) = B(i-1);A(i,i+1) = C(i-1);endM = A\d;syms x;for i = 1:nSx(i) = collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))+M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);digits(4);Sx(i) = vpa(Sx(i));endfor i = 1:ndisp('S(x)=');fprintf(' %s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endS = zeros(1,n);for i = 1:nx = X(i)+;S(i) = Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))+M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3;enddisp('S(i+');disp(' i X(i+ S(i+ ');for i = 1:nfprintf(' %d %.4f %.4f\n',i,X(i)+,S(i));endend输出结果:>> X = [0 1 2 3 4 5 6 7 8 9 10];>> Y = [ ];>> Threch(X,Y,,S(x)=*x - *x^2 - *x^3 + (0,1)S(x)=*x - *x^2 - *x^3 + (1,2)S(x)=*x - *x^2 - *x^3 + (2,3)S(x)=*x^2 - *x - *x^3 + (3,4)S(x)=*x - *x^2 + *x^3 - (4,5)S(x)=*x^2 - *x - *x^3 + (5,6)S(x)=*x - *x^2 + *x^3 - (6,7)S(x)=*x^2 - *x - *x^3 + (7,8)S(x)=*x - *x^2 + *x^3 - (8,9)S(x)=*x - *x^2 + *x^3 - (9,10)S(i+i X(i+ S(i+12345678910[ - *x^3 - *x^2 + *x + , - *x^3 - *x^2 + *x + , - *x^3 - *x^2 + *x + , - *x^3 + *x^2 - *x + , *x^3 - *x^2 + *x - , - *x^3 + *x^2 - *x + , *x^3 - *x^2 + *x - , - *x^3 + *x^2 - *x + , *x^3 - *x^2 + *x - , *x^3 - *x^2 + *x - ]大作业六1、炼钢厂出钢时所用的圣刚睡的钢包,在使用过程中由于钢液及炉渣对包衬耐火材料的侵蚀,使其容积不断增大,经试验,钢包的容积与相应的使用次数的数据如下:(使用次数x,容积y)选用线1=ya对使用最小二乘法数据进行拟合。

相关文档
最新文档