广东工业大学数值计算试卷及答案

广东工业大学试卷用纸，共3页，第2页数值计算引论学院机电工程学院专业机械设计制造及其自动化年级班别2014级（6）班学号3114000271学生姓名刘就杰2016年11月一编写雅可比迭代法求解线性方程组的程序，要求附有算例（20分）。（可能的算例包括基本的验证性算例、方程系数随机生成的一般算例、用于算法对比的比较性算例等，对各算例的结果进行分析。）雅可比迭代法的matlab程序如下functionx=Jacobi(A,b,x0,tol)%雅可比迭代法解线性方程组%A为系数矩阵，b为右端项，x0为初始向量，tol为误差精度sprintf('USAGE:Jacobi(A,b,x0,tol)')D=diag(diag(A));%diag(x)返回由向量x的元素构成的对角矩阵U=triu(A,1);%triu(A)提取矩阵A的上三角部分生成上三角矩阵L=tril(A,-1);%tril(A)提取矩阵A的下三角部分生成下三角矩阵B=-D\(L+U);%B为迭代矩阵dl=D\b;x=B*x0+dl;n=1;whilenorm(x-x0)>=tolx0=x;x=B*x0+dl;n=n+1;endn%n为迭代次数高斯-赛德尔迭代法的matlab程序如下：functionx=Guass_seidel(A,b,x0,tol)%高斯-赛德尔迭代法解线性方程组%A为系数矩阵，b为右端项，x0为初始向量，tol为误差精度sprintf('USAGE:Guass_seidel(A,b,x0,tol)')D=diag(diag(A));%diag(x)返回由向量x的元素构成的对角矩阵U=triu(A,1);%triu(A)提取矩阵A的上三角部分生成上三角矩阵L=tril(A,-1);%tril(A)提取矩阵A的下三角部分生成下三角矩阵G=-(D+L)\U;%G为迭代矩阵dl=(D+L)\b;x=G*x0+dl;n=1;whilenorm(x-x0)>=tolx0=x;x=G*x0+dl;n=n+1;endn%n为迭代次数调用编好的程序求解方程组：A=[5-1-1-1;-110-1-1;-1-15-1;-1-1-110];b=[-4;12;8;34];x0=[0;0;0;0];tol=1e-6;x=Jacobi(A,b,x0,tol)x=Guass_seidel(A,b,x0,tol)实验结果如下：ans=USAGE:Jacobi(A,b,x0,tol)n=20x=1.00002.00003.00004.0000ans=USAGE:Guass_seidel(A,b,x0,t)n=12x=1.00002.00003.00004.0000取相同的初始值达到同样的精度10-6，雅可比迭代需要迭代20次，而高斯-赛德尔迭代法只需12次。实验总结：通过这次实验，对雅可比迭代法以及高斯-赛德尔迭代法求解线性方程组的基本原理有了进一步的理解，同时了解了雅可比和高斯-赛德尔迭代法的优点，即雅可比和高斯-赛德尔在求解线性方程组的过程中具有更快的收敛速度，而高斯-赛德尔比雅可比的收敛速度更快（即取相同的初始值，达到同样精度所需的迭代次数较少）。二编写分段二次拉格朗日插值的程序，要求附有算例（20分）。（对QUOTE在节点0，0.2，0.4，0.6，0.8，1.0上进行插值，求x=0.7处的值，绘出被插值函数与插值函数的图形，予以对比。）建立如下拉格朗日插值函数：functiony=lagrange(x0,y0,x);n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end在matlab中用拉格朗日插值求0.7处的值>>exp(0.7)ans=2.013752707470477>>lagrange(x,y,0.7)ans=2.013751960394443绘出被插值函数与插值函数的图形x=[00.20.40.60.81.0];y=exp(x);x0=[-5:0.001:5];y0=lagrange(x,y,x0);y1=exp(x0);plot(x0,y0,'r')holdonplot(x0,y1,'g')红线为插值函数，绿线是被插值函数，由图像可以知道，在区间（-2，2）是较好拟合的，当超出这个范围后就会偏差越来越大。三编写复化辛普森积分的程序，要求附有算例（20分）。（对定积分QUOTE,计算精度达到QUOTE）复化辛普森积分的程序functionS=bianfuhuasimpson(fx,a,b,eps,M)%变步长复合simpson求积公式%fx--求积函数（函数文件）%a,b--求积区间%eps--计算精度%M--最大允许输出划分数n=1;h=(b-a)/n;T1=h*(feval(fx,a)-feval(fx,b))/2;Hn=h*feval(fx,(a+b)/2);S1=(T1+2*Hn)/3;n=2*n;%最好与倒数第三行保持一致（变步长）whilen<=MT2=(T1+Hn)/2;Hn=0;h=(b-a)/n;forj=1:nx(j)=a+(j-1/2)*h;y(j)=feval(fx,x(j));Hn=Hn+y(j);endHn=h*Hn;S2=(T2+2*Hn)/3;fprintf('n=%2dS2=%-12.9fS2-S1=%-12.9f\n',n,S2,abs(S2-S1));ifabs(S2-S1)<epsbreak;elseT1=T2;S1=S2;n=2*n;endendS=S2;程序执行情况对定积分QUOTE,计算精度达到四编写欧拉法、隐式欧拉法求常微分方程初值问题的程序，要求附有算例（20分）。（对初值问题QUOTEQUOTE采用不同步长（h=0.1,0.01,0.001,0.0001),运行两种算法的程序，并将结果绘制成图形，进行比较、分析。若要解的精度达到QUOTE，应采取什么措施？）程序：%Euler法F='x^2+x-y';a=0;b=0.5;h=0.1;n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1);Y(1)=0;fori=2:n+1x=X(i-1);y=Y(i-1);Y(i)=Y(i-1)+eval(F)*h;end%隐式Euler法Y1=zeros(1,n+1);Y1(1)=0;fori=2:n+1x=X(i);y=Y1(i-1);Y1(i)=(y+x*h^2+h*x)/(h+1);end%准确解temp=[];f=dsolve('Dy=x^2+x-y','y(0)=0','x');df=zeros(1,n+1);fori=1:n+1temp=subs(f,'x',X(i));df(i)=double(vpa(temp));enddisp('步长Euler法隐式Euler法准确值');disp([X',Y',Y1',df']);%画图观察效果figure;plot(X,df,'k-',X,Y,'--r',X,Y1,'.-b');gridon;title('Euler法和隐式Euler法解常微分方程');legend('准确值','Euler法','隐式Euler法');结果为：1.当h=0.1时Euler法和,隐式Euler法与实际值相差较大，误差随x增大而增大。由图可以看出h=0.1时，不能得到精确结果。Euler法精度为，隐式Euler法精度为。2.当h=0.01时Euler法和隐式Euler法与实际值相差不大，误差随x增大而增大，但还是相对准确。由图可以看出h=0.01时，能得到精确结果.Euler法精度为，隐式Euler法精度为。3.当h=0.001时，Euler法和,隐式Euler法与实际值十分相近，误差很小。由图可以看出，三条函数曲线重叠。h=0.001时，能得到精确结果，,.Euler法精度为，隐式Euler法精度为。4.当h=0.0001时Euler法和,隐式Euler法与实际值十分相近，误差很小。由图可以看出，三条函数曲线重叠。h=0.0001时，能得到精确结果，，Euler法精度为，Euler预测-校正法精度为。总体来看，隐式Euler法和Euler法精度没有大的差别，步数每缩小10倍，精度提高10倍。若要解的精度达到QUOTE，可以取h=0.000001，隐式Euler法和Euler法能达到要求。五编写简单迭代法求解非线性方程的程序，要求附有算例（20分）。（对方程QUOTE,求x=0.5附近的根，写出5种迭代格式，至少两种格式收敛；对每种格式的计算结果进行分析、比较，并尽可能把不收敛格式转化为收敛形式）formatlongf=inline('此处用用下面构造高的5种式子替代')disp('x=');x=feval(f,0.5);disp(x);Eps=1E-5;i=1;while1x0=x;i=i+1;x=feval(f,x);disp(x);if~isreal(x)%不是实数不进行迭代disp('出现复数')%提示出现复数break;endifx>1E10;%发散不进行迭代disp('发散')%提示发散break;end;ifabs(x-x0)<Epsbreak;endend;i,x（1）xk+1=(3*xk-xk^2-1)^(1/3)>>Untitled3f=内联函数:f(x)=(3*x-x^2-1)^(1/3)x=0.6299605249474370.7899958937191140.9068993084039620.9648565993178080.9877237574731920.9958404055449220.9986057580992460.9995343879686730.9998446996077400.9999482224823110.9999827396358810.9999942464128830.999998082122915i=13x=0.999998082122915由程序运行结果知道该方法收敛，用了13次就达到求解。（2）xk+1=(3*xk-xk^3-1)^(1/2)>>Untitled3f=内联函数:f(x)=(3*x-x^3-1)^(1/2)x=0.6123724356957940.7794085217018480.9299205949226350.9927793307332860.9999219780951630.9999999908691111.000000000000000i=7x=1.000000000000000由程序运行结果知道，该方法收敛，用了7次就达到求解。相比第一种方式更快，更精确。（3）xk+1=(xk^3+xk^2+1)/3>>Untitled3f=内联函数:f(x)=((x^3+x^2+1)/3)x=0.4583333333333330.4354504243827160.4240619688697000.4186956679577010.4162353170820470.4151217911350620.4146208071262940.4143960160615410.4142952745592360.4142501511564190.4142299447301610.414220897204816i=12x=0.414220897204816由程序运行结果知道该方法