数值分析上机题目4

. . .实验一实验项目：共轭梯度法求解对称正定的线性方程组实验容：用共轭梯度法求解下面方程组(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;while rho10(-11) & k10(-7)&k104 k=k+1;if k=1 p=r;else beta=(r1*r1)/(r*r);p=r1+beta*p;end r=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)实验二1、 实验目的：用复化Simpson方法、自适应复化梯形方法和Romberg方法求数值积分。实验容：计算以下定积分(1) (2) (3) 实验要求：（1）分别用复化Simpson公式、自适应复化梯形公式计算要求绝对误差限为，输出每种方法所需的节点数和积分近似值，对于自适应方法，显示实际计算节点上离散函数值的分布图；（2）分析比较计算结果。程序：syms xf=x6/10-x2+x %定义函数f（x）n=input(输入所求导数阶数:) f2=diff(f,x,n) %求f(x)的n阶导数（1）复化梯形clcclearsyms x%定义自变量xf=inline(x6/10-x2+x,x) %定义函数f(x)=x*exp(x)，换函数时只需换该函数表达式即可f2=inline(3*x4 - 2,x) %定义f(x)的二阶导数,输入程序1里求出的f2即可。f3=-(3*x4 - 2)%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，以便求最大值e=0.5*10(-7) %精度要求值 a=0 %积分下限b=2 %积分上限x1=fminbnd(f3,1,2) %求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值for n=2:1000000 %求等分数n Rn=-(b-a)/12*(b-a)/n)2*f2(x1) %计算余项if abs(Rn)e %用余项进行判断break% 符合要求时完毕endendh=(b-a)/n %求hTn1=0 for k=1:n-1 %求连加和 xk=a+k*h Tn1=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）复化Simpsonclcclearsyms x%定义自变量xf=inline(x6/10-x2+x,x) %定义函数f(x)=x*exp(x)，换函数时只需换该函数表达式即可f2=inline(36*x2,x) %定义f(x)的四阶导数，输入程序1里求出的f2即可f3=-(36*x2)%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，一边求最大值e=5*10(-8) %精度要求值 a=0 %积分下限b=2 %积分上限x1=fminbnd(f3,1,2) %求负的四阶导数的最小值点，也就是求四阶导数的最大值点对应的x值for n=2:1000000 %求等分数n Rn=-(b-a)/180*(b-a)/(2*n)4*f2(x1) %计算余项if abs(Rn)e %用余项进行判断break% 符合要求时完毕endendh=(b-a)/n %求hSn1=0 Sn2=0for k=0:n-1 %求两组连加和 xk=a+k*h xk1=xk+h/2 Sn1=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= 1.6674等分数 n= 24764已知值与计算值的误差 R= -6.3976用Simpson公式计算的结果 Sn= 1.0434等分数 n= 76已知值与计算值的误差 R= -6.0216用复化梯形算法计算的结果 Tn= 0.8985等分数 n= 1119已知值与计算值的误差 R= -6.1665用Simpson公式计算的结果 Sn= 0.9406等分数 n= 8已知值与计算值的误差 R= -6.1245用复化梯形算法计算的结果 Tn= 23.2353等分数 n= 1000000已知值与计算值的误差 R= 16.1704用Simpson公式计算的结果 Sn= 23.2617等分数 n= 10647已知值与计算值的误差 R= 16.1969分析：在处理问题时，复化Simpson要比复化梯度计算速度要快很多。2、实验目的：高斯数值积分方法用于积分方程求解。实验容：线性的积分方程的数值求解，可以被转化为线性代数方程组的求解问题。而线性代数方程组所含未知数的个数，与用来离散积分的数值方法的节点个数相同。在节点数相同的前提下，高斯数值积分方法有较高的代数精度，用它通常会得到较好的结果。对第二类Fredholm积分方程首先将积分区间a,b等分成n份，在每个子区间上离散方程中的积分就得到线性代数方程组。实验要求：分别使用如下方法，离散积分方程中的积分1.复化梯形方法；2.复化辛甫森方法；3.复化高斯方法。求解如下的积分方程，方程的准确解为，并比较各算法的优劣。程序结果：当迭代次数n=1时精确解1.00001.28401.64872.11702.7183复化梯形方法0.8591 1.10321.41651.81882.3354复化辛甫森方法0.9993 1.28321.64762.11562.7165复化高斯方法1.0004 1.2846 1.6495 2.1180 2.7195复化梯形方法的平均误差err=0.247复化辛甫森方法的平均误差err=0.00116复化高斯方法的平均误差err=0.0008当迭代次数n=5时，精确解1.00001.28401.64872.11702.7183复化梯形方法0.99341.27551.63782.10302.7003复化辛甫森方法1.00001.28401.64872.11702.7183复化高斯方法1.00001.28401.64872.11702.7183复化梯形方法的平均误差err=0.00116复化辛甫森方法和复化高斯方法的平均误差err=0可以看出，复化高斯方法得到的结果精度最高，复化辛普森方法比复化梯形方法的精度要高，程序：clear,clca=0;b=1;n=5;figurefun1(a,b,n);hold ona=0;b=1;n=5;figurefun2(a,b,n);hold onfigurefun3(a,b,n);编制m函数：function y=Kfun(t,s)y=2/(exp(1)-1)*exp(t);function y=ffun(t)y=-exp(t);function y=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;for i=1:nx(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;for i=1:nw(2*i)=2/3*h;if in w(2*i+1)=h/3;endendfunction x,w=fhtrapz(a,b,n)h=(b-a)/n;x=a:h:b;for i=1:n+1if i=1|i=n+1 w(i)=h/2;else w(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);V L=eig(CM);x ind=sort(diag(L);V=V(:,ind);w=2*V(:,1).2;x=x*(b-a)/2+(b+a)/2;w=w*(b-a)/2;function fun1(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);for i=1:nB(i)=ffun(x1(i);endfor i=1:nfor j=1:n A(i,j)=w1(j)*Kfun(x1(i),x1(j); endendA=eye(n)-A;y1=(AB);yN=x;for i=1:n1+1yN(i)=ffun(x(i);for j=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(复化值,真实值);function fun2(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);for i=1:nB(i)=ffun(x1(i);endfor i=1:nfor j=1:n A(i,j)=w1(j)*Kfun(x1(i),x1(j); endendA=eye(n)-A;y1=(AB);yN=x;for i=1:n1+1yN(i)=ffun(x(i);for j=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(复化值,真实值);function fun3(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);for i=1:nB(i)=ffun(x1(i);endfor i=1:nfor j=1:n A(i,j)=w1(j)*Kfun(x1(i),x1(j); endendA=eye(n)-A;y1=(AB);yN=x;for i=1:n1+1yN(i)=ffun(x(i);for j=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法提供。程序：子程序function E=Euler(fun,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=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)hold onfunction E=Euler_modify(fun,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=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; for n=1:N x(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,0.1,100)Euler_modify(fun,0,1,0.1,100)x,y=Rk_N4(fun,0,1,0.1,100)plot(x,y,-*r)hold onx1=0:0.1:pi;y1=cos(x1)+sin(x1)plot(x1,y1,-.b)title(误差分析);xlabel(x轴);ylabel(y轴);legend(Euler,Euler改进法,R_K法,精确);axis(0 pi -1.5 1.5);grid on画出图形进行比较：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进行运算，观察计算的结果有什么特点？解的曲线是否有界？解的曲线是不是周期的或趋于某个固定点？（2）在问题允许的围适当改变其中的参数值s, r, b，再选取不同的初始值Y0进行试算，观察并记录计算的结果有什么特点？是否发现什么不同的现象？程序：blzeq函数程序 function ydot=blzeq(t,y)global SIGMA RHO BETA A=-BETA 0 y(2);0 -SIGMA SIGMA;-y(2) RHO -1; ydot=A*y;主程序：clfclcglobal SIGMA RHO BETA SIGMA=10.; RHO=28.; BETA=8/3; axis(0 50 -30 30 -30 30) view(3) hold ontitle(Lorenz Attractor) y0=50,50,50; tfinal=100; t,y=ode23(blzeq,0 tfinal,y0); plot3(y(:,1),y(:,2),y(:,3) （1）实验结果与其分析： 题中的方程与程序中的方程的关系是变量进行了轮换，x换成了y，y换成了z，z换成了x。 原点为原方程的一个奇点，当初始位置稍稍偏离原点如取为0，eps，0，（按原方程中的顺 序，下同）得到的图像如下：分析：这是一个典型的奇怪吸引子的图像，曲线有界，但他不收敛于某一点也不是周期的， 而是在两个位置附近来回的跳跃。（2） 取初始位置分别为20，20，20，60，60，60得到的图像如下： 20，20，20 60，60，60分析：初始变量值相同时，曲线总是被吸引回奇怪吸引子附近作来回跳跃。只变化b取s=10,r=20,b=10初值： 20，20，20 的图像：取s=10,r=20,b=5初值： 20，20，20 的图像：取s=10,r=20,b=8/9初值： 20，20，20 的图像：改变s,r也会有发生类似的情况。15 / 15