数值分析上机题目

--本页仅作为文档封面，使用时请直接删除即可--

--内页可以根据需求调整合适字体及大小本页仅作为文档封面，使用时请直接删除即可--

--内页可以根据需求调整合适字体及大小--数值分析上机题目4(总21页)PAGE实验一实验项目：共轭梯度法求解对称正定的线性方程组实验内容：用共轭梯度法求解下面方程组(1)迭代20次或满足时停止计算。编制程序：储存m文件function[x,k]=CGmethod(A,b)n=length(A);x=2*ones(n,1);r=b-A*x;rho=r'*r;k=0;whilerho>10^(-11)&k<1000k=k+1;ifk==1p=r;elsebeta=rho/rho1;p=r+beta*p;endw=A*p;alpha=rho/(p'*w);x=x+alpha*p;r=r-alpha*w;rho1=rho;rho=r'*r;end运行程序：clear,clcA=[2-100;-13-10;0-14-1;00-15];b=[3-215]';[x,k]=CGmethod(A,b)运行结果：x=(2),A是1000阶的Hilbert矩阵或如下的三对角矩阵，A[i,i]=4，A[i,i-1]=A[i-1,i]=-1，i=2,3,..,nb[1]=3,b[n]=3,b[i]=2，i=2,3,…,n-1迭代10000次或满足时停止计算。编制程序：储存m文件function[x,k]=CGmethod_1(A,b)n=length(A);x(1:n,1)=0;r=b-A*x;r1=r;k=0;whilenorm(r1,1)>10^(-7)&k<10^4k=k+1;ifk==1p=r;elsebeta=(r1'*r1)/(r'*r);p=r1+beta*p;endr=r1;w=A*p;alpha=(r'*r)/(p'*w);x=x+alpha*p;r1=r-alpha*w;end运行程序：clear,clcn=1000;A=hilb(n);b=sum(A')';[x,k]=CGmethod(A,b)实验二实验目的：用复化Simpson方法、自适应复化梯形方法和Romberg方法求数值积分。实验内容：计算下列定积分(1)(2)(3)实验要求：（1）分别用复化Simpson公式、自适应复化梯形公式计算要求绝对误差限为，输出每种方法所需的节点数和积分近似值，对于自适应方法，显示实际计算节点上离散函数值的分布图；（2）分析比较计算结果。程序：symsxf=x^6/10-x^2+x%定义函数f（x）n=input('输入所求导数阶数:')f2=diff(f,x,n)%求f(x)的n阶导数（1）复化梯形clcclearsymsx%定义自变量xf=inline('x^6/10-x^2+x','x')%定义函数f(x)=x*exp(x)，换函数时只需换该函数表达式即可f2=inline('3*x^4-2','x')%定义f(x)的二阶导数,输入程序1里求出的f2即可。f3='-(3*x^4-2)'%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，以便求最大值e=*10^(-7)%精度要求值a=0%积分下限b=2%积分上限x1=fminbnd(f3,1,2)%求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值forn=2:1000000%求等分数nRn=-(b-a)/12*((b-a)/n)^2*f2(x1)%计算余项ifabs(Rn)<e%用余项进行判断break%符合要求时结束endendh=(b-a)/n%求hTn1=0fork=1:n-1%求连加和xk=a+k*hTn1=Tn1+f(xk)endTn=h/2*((f(a)+2*Tn1+f(b)))z=exp(2)R=Tn-z%求已知值与计算值的差stem(xk,Tn1);fprintf('用复化梯形算法计算的结果Tn=')disp(Tn)fprintf('等分数n=')disp(n)%输出等分数fprintf('已知值与计算值的误差R=')disp(R)（2）复化Simpsonclcclearsymsx%定义自变量xf=inline('x^6/10-x^2+x','x')%定义函数f(x)=x*exp(x)，换函数时只需换该函数表达式即可f2=inline('36*x^2','x')%定义f(x)的四阶导数，输入程序1里求出的f2即可f3='-(36*x^2)'%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，一边求最大值e=5*10^(-8)%精度要求值a=0%积分下限b=2%积分上限x1=fminbnd(f3,1,2)%求负的四阶导数的最小值点，也就是求四阶导数的最大值点对应的x值forn=2:1000000%求等分数nRn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1)%计算余项ifabs(Rn)<e%用余项进行判断break%符合要求时结束endendh=(b-a)/n%求hSn1=0Sn2=0fork=0:n-1%求两组连加和xk=a+k*hxk1=xk+h/2Sn1=Sn1+f(xk1)Sn2=Sn2+f(xk)endSn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a))+f(b))%因Sn2多加了k=0时的值，故减去f(a)z=exp(2)R=Sn-z%求已知值与计算值的差fprintf('用Simpson公式计算的结果Sn=')disp(Sn)fprintf('等分数n=')disp(n)fprintf('已知值与计算值的误差R=')disp(R)结果：用复化梯形算法计算的结果Tn=等分数n=24764已知值与计算值的误差R=用Simpson公式计算的结果Sn=等分数n=76已知值与计算值的误差R=用复化梯形算法计算的结果Tn=等分数n=1119已知值与计算值的误差R=用Simpson公式计算的结果Sn=等分数n=8已知值与计算值的误差R=用复化梯形算法计算的结果Tn=等分数n=1000000已知值与计算值的误差R=用Simpson公式计算的结果Sn=等分数n=10647已知值与计算值的误差R=分析：在处理问题时，复化Simpson要比复化梯度计算速度要快很多。2、实验目的：高斯数值积分方法用于积分方程求解。实验内容：线性的积分方程的数值求解，可以被转化为线性代数方程组的求解问题。而线性代数方程组所含未知数的个数，与用来离散积分的数值方法的节点个数相同。在节点数相同的前提下，高斯数值积分方法有较高的代数精度，用它通常会得到较好的结果。对第二类Fredholm积分方程首先将积分区间[a,b]等分成n份，在每个子区间上离散方程中的积分就得到线性代数方程组。实验要求：分别使用如下方法，离散积分方程中的积分1.复化梯形方法；2.复化辛甫森方法；3.复化高斯方法。求解如下的积分方程，方程的准确解为，并比较各算法的优劣。程序结果：当迭代次数n=1时精确解复化梯形方法复化辛甫森方法复化高斯方法复化梯形方法的平均误差err=复化辛甫森方法的平均误差err=复化高斯方法的平均误差err=当迭代次数n=5时，精确解复化梯形方法复化辛甫森方法复化高斯方法复化梯形方法的平均误差err=复化辛甫森方法和复化高斯方法的平均误差err=0可以看出，复化高斯方法得到的结果精度最高，复化辛普森方法比复化梯形方法的精度要高，程序：clear,clca=0;b=1;n=5;figurefun1(a,b,n);holdona=0;b=1;n=5;figurefun2(a,b,n);holdonfigurefun3(a,b,n);编制m函数：functiony=Kfun(t,s)y=2/(exp(1)-1)*exp(t);functiony=ffun(t)y=-exp(t);functiony=Fexc(t)%精确解y=exp(t);function[x,w]=fhgauss(a,b,n)h=(b-a)/n;x1=a:h:b;x=zeros(1,2*n);w=x;fori=1:n[x(2*i-1:2*i),w(2*i-1:2*i)]=GaussLegendre(x1(i),x1(i+1),2);endfunction[x,w]=fhsimpson(a,b,n)h=(b-a)/n;x=a:h/2:b;w=x;w(1)=h/6;w(2*n+1)=h/6;fori=1:nw(2*i)=2/3*h;ifi<nw(2*i+1)=h/3;endendfunction[x,w]=fhtrapz(a,b,n)h=(b-a)/n;x=a:h:b;fori=1:n+1ifi==1||i==n+1w(i)=h/2;elsew(i)=h;endendfunction[x,w]=GaussLegendre(a,b,n)i=1:n-1;c=i./sqrt(4*i.^2-1);CM=diag(c,1)+diag(c,-1);[VL]=eig(CM);[xind]=sort(diag(L));V=V(:,ind)';w=2*V(:,1).^2;x=x*(b-a)/2+(b+a)/2;w=w*(b-a)/2;functionfun1(a,b,n)[x1,w1]=fhtrapz(a,b,n);n1=4;n=n+1;h=(b-a)/n1;x=a:h:b;y0=Fexc(x);A=zeros(n,n);B=zeros(n,1);fori=1:nB(i)=ffun(x1(i));endfori=1:nforj=1:nA(i,j)=w1(j)*Kfun(x1(i),x1(j));endendA=eye(n)-A;y1=(A\B)';yN=x;fori=1:n1+1yN(i)=ffun(x(i));forj=1:nyN(i)=yN(i)+w1(j)*Kfun(x(i),x1(j))*y1(j);endendfprintf('数值解:\n')yNfprintf('精确解:\n')y0plot(x,yN,'x',x,y0,'.');h=legend('复化值','真实值');functionfun2(a,b,n)[x1,w1]=fhsimpson(a,b,n);n=2*n+1;n1=50;h=(b-a)/n1;x=a:h:b;y0=Fexc(x);A=zeros(n,n);B=zeros(n,1);fori=1:nB(i)=ffun(x1(i));endfori=1:nforj=1:nA(i,j)=w1(j)*Kfun(x1(i),x1(j));endendA=eye(n)-A;y1=(A\B)';yN=x;fori=1:n1+1yN(i)=ffun(x(i));forj=1:nyN(i)=yN(i)+w1(j)*Kfun(x(i),x1(j))*y1(j);endendfprintf('数值解:\n')yNfprintf('精确解:\n')y0plot(x,yN,'x',x,y0,'.');h=legend('复化值','真实值');functionfun3(a,b,n)[x1,w1]=fhgauss(a,b,n);n=2*n;n1=4;h=(b-a)/n1;x=a:h:b;y0=Fexc(x);A=zeros(n,n);B=zeros(n,1);fori=1:nB(i)=ffun(x1(i));endfori=1:nforj=1:nA(i,j)=w1(j)*Kfun(x1(i),x1(j));endendA=eye(n)-A;y1=(A\B)';yN=x;fori=1:n1+1yN(i)=ffun(x(i));forj=1:nyN(i)=yN(i)+w1(j)*Kfun(x(i),x1(j))*y1(j);endendfprintf('数值解:\n')yNfprintf('精确解:\n')y0plot(x,yN,'x',x,y0,'.');h=legend('复化值','真实值');图一图二图三实验三1、对常微分方程初值问题分别使用Euler显示方法、改进的Euler方法和经典RK法和四阶Adams预测-校正算法，求解常微分方程数值解，并与其精确解进行作图比较。其中多步法需要的初值由经典RK法提供。程序：子程序functionE=Euler(fun,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;forn=1:Nx(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(fun,x(n),y(n));endx1=xy1=yplot(x1,y1,'k')holdonfunctionE=Euler_modify(fun,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;forn=1:Nx(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n));y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n))+feval(fun,x(n+1),z0));endx2=xy2=yplot(x2,y2,'g')function[x,y]=Rk_N4(f,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;forn=1:Nx(n+1)=x(n)+h;k1=h*feval(f,x(n),y(n));k2=h*feval(f,x(n)+1/2*h,y(n)+1/2*k1);k3=h*feval(f,x(n)+1/2*h,y(n)+1/2*k2);k4=h*feval(f,x(n)+h,y(n)+k3);y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);end运行以下程序clcclearEuler('fun',0,1,,100)Euler_modify('fun',0,1,,100)[x,y]=Rk_N4('fun',0,1,,100)plot(x,y,'-*r')holdonx1=0::pi;y1=cos(x1)+sin(x1)plot(x1,y1,'')title('误差分析');xlabel('x轴');ylabel('y轴');legend('Euler','Euler改进法','R_K法','精确');axis([0pi]);gridon画出图形进行比较：2、实验目的：Lorenz问题与混沌实验内容：考虑著名的Lorenz方程其中s,r,b为变化区域有一定限制的实参数。该方程形式简单，表面上看并无惊人之处，但由该方程揭示出的许多现象，促使“混沌”成为数学研究的崭新领域，在实际应用中也产生了巨大的影响。实验方法：先取定初值Y0=(x,y,z)=(0,0,0)，参数s=10,r=28,b=8/3，用MATLAB编程数值求解，并与MATLAB函数ods45的计算结果进行对比。实验要求：（1）对目前取定的参数值s,r和b，选取不同的初值Y0进行运算，观察计算的结果有什么特点解的曲线是否有界解的曲线是不是周期的或趋于某个固定点