山建大数值分析实验报告Word格式文档下载.docx

1、n-1 for i=k+1:n a（i,k）=a（i,k）/a（k,k）; b（i）=b（i）-a（i,k）*b（k）; for j=k+1: if a（k,k）=0主元素为零，消去法无法继续） ; break; else a（i,j）=a（i,j）-a（i,k）*a（k,j）; endb（n）=b（n）/a（n,n）;for i=（n-1）:-1:1 w=0; for j=（i+1）: w=w+a（i,j）*b（j）; b（i）=（b（i）-w）/a（i,i）; y=b;（b）高斯列主元消去法 function z=gauss2（a,b,ep） %高斯列主元素消元法if nargin=2 e

2、p=0.000001 p=a（k,k）;I=k; for i=k: if abs（a（i,k）abs（p） p=a（i,k）;I=i; if p 2 -1 0 0 -1 -3 -2 0 -1 3 -2 0 0 0 -3 5 b = 6 1 0 gauss1（a,b）方程组的解为ans = 35/12 -1/6 -41/24 -33/40 （6）实验体会：主元消去法和高斯消去法的确是两个非常锻炼人编程的方法，在编写程序时，需要使用的大量的循环和分支结构，但无论是高斯消去法还是高斯列主元法，它们的原理还算不难理解，通过变成能够较好的理解它们。实验二 解线性方程组的迭代法2.1实验目的 掌握解线性方

3、程组的雅可比迭代和高斯-塞德尔迭代算法； 培养编程与上机调试能力.2.2实验要求： （1）选择一种计算机语言（Matlab）设计出雅可比（Jacobi）Gauss-Seidel、SOR迭代法，迭代法的算法程序，并且选择不同的迭代次数，观察输出结果； （2）利用Matlab求方程组2.3 实验内容计算书上的习题P61例2.5.1a=20 2 3;1 8 1;2 -3 15; b=24;12;30; x0=0;0;0;（2）对应程序：Jacobi迭代法：function X=jacobi（a,b,X0,ep）%Jacobi迭代法求解方程组if nargin=3 ep=1.0e-6;elseif n

4、argin=ep）&（k=1000） X0=X; X=B*X0+f; k=k+1;disp（迭代次数为kreturn高斯-塞德尔（Gauss-Seidel）：function z=gauss2（a,b,ep）gauss2（a,b,ep）（3）实验结果：迭代次数为k = 6 0. 33 1. 33 2. 67高斯消元法： 0. 15 1. 94 2. 37（4）实验体会： Jacobi迭代法和高斯消元法都能很好的解决方程组的求解问题，在上机程中遇到了也遇到了不少问题，但最后在老师的悉心辅导下都得到了很好的解答，这两个程序使我明白了要变出好的程序就需要我们积极思考问题，勇于发现和解决问题。实验三

5、矩阵特征值问题计算3.1实验目的 掌握求矩阵的特征值和主特征向量的幂法；3.2实验要求（1） 选择一种计算机语言设计出幂法求主特征值和相应特征向量的程序，并且选择不同的初值，观察所需的迭代次数和迭代结果.（2） 利用Matlab求特征值和特征向量 调用格式1： eig（A） %得到特征值列向量调用格式2：，其中为由特征列向量构成的方阵，为由特征值构成的对角阵. %得到特征值和所对应的特征向量3.3 实验内容P81例3.1.1A=2 4 6;3 9 15;4 16 36 x0=1;1;1 乘幂法function y=chm（a,x0,k）%乘幂法求主特征值及特征向量if nargin=eps）&

6、10000） a=G; R=R1;for i=2: if abs（a（i,j）w w=abs（a（i,j）; p=i; q=j; if a（p,p）=a（q,q） theta=sign（a（p,q）*pi/4; G（p,p）=a（p,p）*（cos（theta）2+a（q,q）*（sin（theta）2+a（p,q）*sin（2*theta）; G（q,q）=a（p,p）*（sin（theta）2+a（q,q）*（cos（theta）2-a（p,q）*sin（2*theta）; G（p,q）=（a（q,q）-a（p,p）*（cos（theta）*（sin（theta）+a（p,q）*（cos（

7、2*theta）; G（q,p）=G（p,q）; for i=1: R1（i,p）=R（i,p）*cos（theta）+R（i,q）*sin（theta）; R1（i,q）=-R（i,p）*sin（theta）+R（i,q）*cos（theta）; if （i=p）&（i=q） G（p,i）=a（i,p）*（cos（theta）+a（i,q）*（sin（theta）; G（i,p）=G（p,i）; G（q,i）=-a（i,p）*（sin（theta）+a（i,q）*（cos（theta）; G（i,q）=G（q,i）; else gasi=（a（p,p）-a（q,q）/（2*a（p,q）; t

8、=sign（gasi）*（-abs（gasi）+sqrt（1+gasi2）; G（p,p）=a（p,p）*（1/（1+t2）+a（q,q）*（t2/（1+t2）+a（p,q）*2*（t/（1+t2）; G（q,q）=a（p,p）*（t2/（1+t2）+a（q,q）*（1/（1+t2）-a（p,q）*2*（t/（1+t2）; G（p,q）=（a（q,q）-a（p,p）*（t/（1+t2）+a（p,q）*（1/（1+t2）-（t2/（1+t2）; R1（i,p）=R（i,p）*（1/sqrt（1+t2）+R（i,q）*（t/sqrt（1+t2）; R1（i,q）=-R（i,p）*（t/sqrt（1

9、+t2）+R（i,q）*（1/sqrt（1+t2）; G（p,i）=a（i,p）*（1/sqrt（1+t2）+a（i,q）*（t/sqrt（1+t2）; G（q,i）=-a（i,p）*（t/sqrt（1+t2）+a（i,q）*（1/sqrt（1+t2）; SS=SS+G（i,j）2;G的对角线即为近似的特征值 GR的列向量为相应的近似特征向量）, RJeig（A,10e-9）乘幂法：m = 43. 56j = 3lanmda =u1 = 0. 50 0. 10 1. 00经典Jacobi法：G的对角线即为近似的特征值G = 0. 12 0. 20 0. 20 0. 20 2. 88 -0. 7

10、5 0. 20 -0. 75 44. 99R的列向量为相应的近似特征向量R = 0. 54 0. 01 0. 31 -0. 62 0. 21 0. 06 0. 92 -0. 39 0. 44 本次程序是由老师提供的，但是同学们都认真阅读过了，我发现此程序是相当复杂的，要对矩阵进行迭代等很多操作，我还了解到我们在处理较大的问题是必须要使用到程序，因此让我对程序产生的浓厚的兴趣，同时我们也认识到我们所编写的程序稳定性很差，因此我们还需要更多的练习。实验四 插值法4.1 实验目的1 掌握插值法的基本思路和步骤；2 培养编程与上机调试能力。4.2 实验要求用Matlab和插值中的某种具体算法编写代码并

11、执行，完成解决具体问题。4.3 实验内容Matlab SplineToolsx = 0 0.2000 0.4000 0.6000 0.8000y = 1.0000 1.2214 1.4918 1.8221 2.2255x0 = 0.1500function s=Newton（x,y,x0,nn）%nn为Newton插值多项式的次数nx=length（x）;ny=length（y）;if nx=ny warning （向量x与y的长度应该相同 return;m=length（x0）;m yy=y; t=0; kk=1; while（kk=nn） kk=kk+1; for k=kk:nx yy（k

12、）=（yy（k）-yy（kk-1）/（x（k）-x（kk-1）; t=yy（1）; for k=2:nn+1 u=1.0; j=1; while（j=ep area=0.0; n=n+1; h=h/2; area=area+feval（fun,h*（2*i-1）+a）; t（n+1,1）=0.5*t（n,1）+area*h; m=2*m; if n4 n-j t（i,j+1）=（4（j）*t（i+1,j）-t（i,j）/（4（j）-1）; t1=t（n-4,4）; t2=t（n-3,4）;end 用Romberg序列求得积分近似值为s=t2;function y=fun（x）if x=0 y=

13、1;y=sin（x）/xromberg1（fune-6） 0. 26 以前不知道用程序怎么求积分，但通过本次实验我知道了利用程序是可以求积分的，它主要用的的有限次的迭代，将连续的无限多点的问题有线化，然后就可以利用计算机来解决积分问题，本程序由于比较复杂，所以老师给出了，虽然程序是由老师给出的，但是我们在调用该程序时也基本明白了程序实现的，在实验过程中我们仔细阅读了程序，所以基本上可以再现，但我个人觉得这还远远不够，我们需要自己独立思考学习编写，不断的提高自己的编程能力。实验七常微分方程数值解法7.1实验目的 掌握常微分方程数值解的常用算法；求解常微分方程初值、边值问题的解.7.2实验题目：下

14、述方面的相应习题（1）用改进的欧拉公式，求解常微分方程初值问题的解 （2）用四阶龙格-库塔公式解初值问题（3）求解常微分方程边值问题的解7.3实验要求（1） 选择一种计算机语言设计出改进欧拉法和四阶龙格-库塔法方法求解常微分方程初值问题的程序，观察运行结果.（2） 利用Matlab求解常微分方程初值问题 函数dsolve（）用于求解微分方程.Dy表示：dy/dt（t 为缺省的自变量），Dny表示y对t的n阶导数. Matlab6.5环境下操作如下： y=dsolve（Dy=y*y,y（0）=1） %求解题目1Dy=y/ty（2.0）=1） %求解题目2 （3）ode23， ode45（1） 利

15、用最小二乘法拟合通过改进欧拉法求出微分方程的一系列数值解的近似函数方程.并利用Matlab的绘图功能画出函数的曲线。（1）题目1：y=x2+y2,0x1, y（0）=0, h=0.1题目2： y=3*y/（1+x）,01, y（0）=1, h=0.2（4）对应程序1：function s=Eulerc（fun,x0,xn,y0,h）%h为步长，改进的求微分方程数值解的Euler公式。5n=（xn-x0）/h; t1=y0+h*feval（fun,x0,y0）; x0=x0+h; t2=y0+h*feval（fun,x0,t1）; y0=（t1+t2）/2; y（i）=y0;s=y;对应程序2：function s=RK（fun,x0,xn,y0,h）%h为步长n=（x