数值分析课程设计大作业

课程设计2013年07月20日设计题目《数值分析》课程设计学生姓名****学号####专业班级指导教师设计题目《数值分析》课程设计成绩课程设计主要内容实验一：三题全做；实验二：四题全做；实验三：3.13.2已做；实验四4.14.24.34.4已做；实验五：5.15.25.35.4已做；实验六：6.26.36.4已做；实验七：7.37.57.67.8已做；实验八：8.28.8已做。共计26道题。部分题通过matlab中notebook实现程序，notebook使用非常方便。由于篇幅有限部分题调用的函数源代码没有给出，本作业所使用的m文件均在附录中指导教师评语建议：从学生的工作态度、工作量、设计（论文）的创造性、学术性、实用性及书面表达能力等方面给出评价。签名：20年月日1.1水手、猴子和椰子问题算法分析:设椰子起初的数目为0p,第一至第五次猴子在夜里藏椰子后，椰子的数目分别为0p，1p，2p，3p，4p，再设最后每个人分得x个椰子，由题意得：15541(1),0,1,2,3,4.(1)=5155kkppkxppx所以利用逆向递推方法求解：n=input('n=');forx=1:np=5*x+1;fork=1:5p=5*p/4+1;endifp==fix(p)breakendenddisp([x,p])执行代码后得：n=102315621（输入n=1000000）即最后每个人分得1023个椰子，椰子总数为156211.2当0,1,2,,100n时，选择稳定的算法计算积分10d10nxnxeIxe由1100(1)1110010101,1010110(1)10xxnxnxnxnnnxeIIdxeeeIIdxedxeen得0111(1)1011[(1)],100,99,...,1.10nnnIIIeInn由上式可知求nI时，1nI的误差的影响被缩小了。n=100时100I的近似值为0。matlab代码为fprintf('稳定算法:\n')y0=0;n=100;plot(n,y0,'r*');holdonfprintf('y[100]=%10.6f',y0);while(1)y1=1/10*[(1-exp(-n))/n-y0];fprintf('y[%10.0f]=%10.6f',n-1,y1);plot(n-1,y1,'r*')if(n=1)break;endy0=y1;n=n-1;ifmod(n,3)==0,fprintf('\n'),end,end(具体值已省略)编程实现得下图。010203040506070809010000.010.020.030.040.050.06由图可知，该算法是稳定的。1.3绘制静态和动态的Koch分形曲线Koch曲线程序koch.mfunctionkoch(a1,b1,a2,b2,n)%koch(0,0,9,0,3)%a1,b1,a2,b2为初始线段两端点坐标,n为迭代次数%例如a1=0;b1=0;a2=9;b2=0;n=3;%第i-1次迭代时由各条线段产生的新四条线段的五点横、纵坐标存储在数组A、B中[A,B]=sub_koch1(a1,b1,a2,b2);fori=1:nforj=1:length(A)/5;w=sub_koch2(A(1+5*(j-1):5*j),B(1+5*(j-1):5*j));fork=1:4[AA(5*4*(j-1)+5*(k-1)+1:5*4*(j-1)+5*(k-1)+5),BB(5*4*(j-1)+5*(k-1)+1:5*4*(j-1)+5*(k-1)+5)]=sub_koch1(w(k,1),w(k,2),w(k,3),w(k,4));endendA=AA;B=BB;endplot(A,B)holdonaxisequal%由以(ax,ay),(bx,by)为端点的线段生成新的中间三点坐标并把这五点横、纵坐标依次分别存%储在数组A,B中function[A,B]=sub_koch1(ax,ay,bx,by)cx=ax+(bx-ax)/3;cy=ay+(by-ay)/3;ex=bx-(bx-ax)/3;ey=by-(by-ay)/3;L=sqrt((ex-cx).^2+(ey-cy).^2);alpha=atan((ey-cy)./(ex-cx));if(ex-cx)0alpha=alpha+pi;enddx=cx+cos(alpha+pi/3)*L;dy=cy+sin(alpha+pi/3)*L;A=[ax,cx,dx,ex,bx];B=[ay,cy,dy,ey,by];%把由函数sub_koch1生成的五点横、纵坐标A,B顺次划分为四组,分别对应四条折线段中%每条线段两端点的坐标,并依次分别存储在4*4阶矩阵k中,k中第i(i=1,2,3,4)行数字代表第%i条线段两端点的坐标functionw=sub_koch2(A,B)a11=A(1);b11=B(1);a12=A(2);b12=B(2);a21=A(2);b21=B(2);a22=A(3);b22=B(3);a31=A(3);b31=B(3);a32=A(4);b32=B(4);a41=A(4);b41=B(4);a42=A(5);b42=B(5);w=[a11,b11,a12,b12;a21,b21,a22,b22;a31,b31,a32,b32;a41,b41,a42,b42];%调用函数得到下图n=5;i=0;whileinfigure(i+1);koch(0,0,3,0,i)i=i+1;pause(1)end将pause(1)去掉可得静态图00.511.522.53-0.500.511.500.511.522.53-0.500.511.500.511.522.53-0.500.511.500.511.522.53-0.500.511.500.511.522.53-0.500.511.52.1小行星轨道问题为了确定方程中的5个待定系数，需要将上述5个点的坐标代入上面的方程22123452221axaxyayaxay,得：22112113141512212222324252221323333435322142443444542215255354555ax2axyay2ax2ay1ax2axyay2ax2ay1ax2axyay2ax2ay1ax2axyay2ax2ay1ax2axyay2ax2ay1将这一包含5个未知数的线性方程组，写成矩阵的形式2211111122222222223333332244444422555555x2xyy2x2yx2xyy2x2yx2xyy2x2yx2xyy2x2yx2xyy2x2y12345aaaaa=11111AXb（1）求解这一线性方程组，即可得到曲线方程的系数X0=[5360558460628596666268894];Y0=[606211179169542349268894];A=zeros(5);X0(1);fori=1:5A(i,1)=X0(i)^2;A(i,2)=2*X0(i)*Y0(i);A(i,3)=Y0(i)^2;A(i,4)=2*X0(i);A(i,5)=2*Y0(i);end;formatlongg;AA=2873496025649907020367478441072101212434175716001307048680124970041116920223539512538812131422972287438116125718339084443822244313204740855187406413332446984474638323694927664724746383236137788137788B=[-1-1-1-1-1]';formatlongg;x=A\Bx=-8.06820280371841e-011-7.63620099622306e-011-3.0801152978055e-010-8.89025159419867e-0062.02829368401655e-005(2)用Lu分解法解可得formatlongg;A=[28734960256499070203674784410721012124;3417571600130704868012497004111692022358;3951253881213142297228743811612571833908;4443822244313204740855187406413332446984;474638323694927664724746383236137788137788];B=[-1-1-1-1-1]';[L,U,flag]=LU_Decom(A),formatlongg;x=U\(L\B)flag=OKx=-8.06820280370254e-011-7.63620099622621e-011-3.08011529780421e-010-8.89025159420231e-0062.02829368401615e-005jacobi迭代法：jacobic(A)因为谱半径不小于1，所以Jacobi迭代序列发散，谱半径SRH和B的所有特征值H如下：SRH=4.41963931714337ans=-4.419639317143371.54532928016960.7571387636487321.118388956505330.9987823168197GSC(A)因为谱半径不小于1，所以Gauss-Seidel迭代序列发散，谱半径SRH和B的所有特征值H如下：SRH=1.12218280703645ans=00.09140470453603941.10703104191813+0.183783907450762i1.10703104191813-0.183783907450762i0.997482236058482.2(1)用Gauss列主元消去法、Gauss按比例列主元消去法、Cholesky分解求解下列线性方程组，并彼此互相验证。(2)判断用Jacobi迭代法、Gauss-Seidel迭代法、SOR法(分别取0.8,1.2,1.3,1.6)解下列线性方程组的收敛性.若收敛，再用Jacobi迭代法、Gauss-Seidel迭代法、SOR法(分别取0.8,1.2,1.3,1.6)分别解线性方程组，并比较各种方法的收敛速度.1234124134123421,333,2965,36197.xxxxxxxxxxxxxx(3)用Cholesky分解求解下列线性方程组12345678424024000221213206411418356200216143323218122410394334411142202531011421500633421945xxxxxxxx（1）A=[1-121;-130-3;209-6;1-3-619];b=[1357]';[RA,RB,n,x]=liezy(A,b),[RA,RB,n,x]=bilizy(A,b),cholesky(A,b)列主元因为RA=RB=n,所以此方程组有唯一解.RA=4RB=4n=4x=-8.0000000000000050.3333333333333323.6666666666666682.000000000000000比例主元因为RA=RB=n,所以此方程组有唯一解.