山建大数值分析实验报告

1、山东建筑大学数值计算A实验报告r=i二级学院：理学院专 业：信息与计算科学指导教师：_XXX班级学号：信计XX 20081212XX姓 名：XXXX实验一高斯消去法题目：用选主元素法和高斯消去法求解下列方程组：已知方程增广矩阵： a = TOC o 1-5 h z 2-1006-1-3-201-13-20000-351原理：高斯消去法：相对于约当消去法而言，总的来说就是将增广矩阵化为下三角阵。顺序选主元素法:相对于高斯消去法的唯一不同是要先按当前要排列元素所在列 大小进行排列。设计思想：高斯消去法：首先让每一行的元素除以该行的主对角线元素。然后利用此行使位 于下一行主对角线以前的元素变为0，依

2、次类推。选主元素法：在高斯消去法的基础上，每次进行化上三角阵之前，重新排列各方 程的位置。对应程序：(a)高斯消去法function y=gauss1(a,b)%高斯顺序消去法m,n=size(a);if m=ndisp(输入错误，系数矩证阵只能是方阵)endif n=length(b)disp(输入错误,常数项的个数应与方程的个数相同)endfor k=1:n-1for i=k+1:na(i,k)=a(i,k)/a(k,k);b(i)=b(i)-a(i,k)*b(k);for j=k+1:nif a(k,k)=0disp(主元素为零，消去法无法继续);break;elsea(i,j)=a(i

3、,j)-a(i,k)*a(k,j);endendendendb(n)=b(n)/a(n,n);for i=(n-1)：-1:1w=0;for j=(i+1):nw=w+a(i,j)*b(j);endb(i) = (b(i)-w)/a(i,i);endy=b;(b)高斯列主元消去法function z=gauss2(a,b,ep)%高斯列主元素消元法if nargin=2ep=0.000001endm,n=size(a);if m=ndisp(输入错误，系数矩证阵只能是方阵)endif n=length(b)disp(输入错误,常数项的个数应与方程的个数相同) endfor k=1:n-1p=a

4、(k,k);I=k;for i=k:nif abs(a(i,k)abs(p)p=a(i,k);I=i;endendif pa =2-100-1-3-20-13-2000-35 b =6101 gauss1(a,b) 方程组的解为ans =35/12-1/6-41/24-33/40实验体会：主元消去法和高斯消去法的确是两个非常锻炼人编程的方法，在编写程序时，需要使用 的大量的循环和分支结构，但无论是高斯消去法还是高斯列主元法，它们的原理还算不难理 解，通过变成能够较好的理解它们。实验二解线性方程组的迭代法2.1实验目的掌握解线性方程组的雅可比迭代和高斯-塞德尔迭代算法；培养编程与上机调试能力.2

5、.2实验要求：选择一种计算机语言(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;对应程序：Jacobi迭代法：function X=jacobi(a,b,X0,ep)%Jacobi迭代法求解方程组if nargin=3ep=1.0e-6;elseif nargin=ep)&(kabs(p)p=a(i,k);I=i;endendif

6、p=epz=0;endif I=kfor j=k:nw=a(k,j); a(k,j)=a(I,j);a(I,j)=w;endu=b(k);b(k)=b(I); b(I)=u;endfor i=k+1:na(i,k)=a(i,k)/a(k,k);b(i)=b(i)-a(i,k)*b(k);for j=k+1:na(i,j)=a(i,j)-a(i,k)*a(k,j);endendendb(n)=b(n)/a(n,n);for i=(n-1)：-1:1w=0;for j=(i+1):nw=w+a(i,j)*b(j);endb(i) = (b(i)-w)/a(i,i);end% disp(方程组的解为

7、)；z=b;gauss2(a,b,ep)实验结果：Jacobi迭代法：迭代次数为k =6ans =0.767363020833331.138413958333332.12535579166667高斯消元法：ans =0.767353807320151.138409760201942.12536811106437实验体会：Jacobi迭代法和高斯消元法都能很好的解决方程组的求解问题，在上机程中遇到了也 遇到了不少问题，但最后在老师的悉心辅导下都得到了很好的解答，这两个程序使我明白了 要变出好的程序就需要我们积极思考问题，勇于发现和解决问题。实验三矩阵特征值问题计算3.1实验目的掌握求矩阵的特征值

8、和主特征向量的幕法；培养编程与上机调试能力.3.2实验要求选择一种计算机语言设计出幕法求主特征值和相应特征向量的程序，并且选择不同 的初值，观察所需的迭代次数和迭代结果.利用Matlab求特征值和特征向量调用格式1： eig(A)%得到特征值列向量调用格式2： D, X eig(A)，其中D为由特征列向量构成的方阵，X为由特征值构 成的对角阵.得到特征值和所对应的特征向量3.3实验内容计算书上的习题题目：P81 例 3.1.1A=2 4 6;3 9 15;4 16 36x0=1;1;1对应程序：乘幕法function y=chm(a,x0,k)%乘幕法求主特征值及特征向量if nargin=e

9、ps)&(kww=abs(a(i,j);p=i;q=j;endendendif 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)”2a(p,q)*sin(2*theta); G(p,q) = (a(q,q)-a(p,p)*(cos(theta)*(sin(theta)+a(p,q)*(cos(2*theta);G(q,p)=G(p,q);for

10、 i=1:nR1(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);endendelsegasi=(a(p,p)-a(q,q)/(2*a(p,q);t=sign(gasi)*(-abs(gasi)+sqrt(1+

11、gas2);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);G(q,p)=G(p,q);for i=1:nR1(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+t2)+R(i,q)*(

12、1/sqrt(1+t2);if (i=p)&(i=q)G(p,i)=a(i,p)*(1/sqrt(1+t2)+a(i,q)*(t/sqrt(1+t2);G(i,p)=G(p,i);G(q,i)=-a(i,p)*(t/sqrt(1+t2)+a(i,q)*(1/sqrt(1+t2);G(i,q)=G(q,i);endendendSS=0;for i=2:nfor j=1:i-1SS=SS+G(i,j)2;endendk=k+1;enddisp（G的对角线即为近似的特征值）Gdisp（R的列向量为相应的近似特征向量），RJeig（A,10e-9）（3）实验结果：乘幕法：m =43.879997817

13、19956j =3lanmda =43.87999781719956u1 =0.185867901919500.446032888095101.00000000000000经典Jacobi法：G的对角线即为近似的特征值G =0.038388746529120.000000000000200.000000000137200.000000000000202.33485967315988-0.000000863196750.00000000013720 -0.00000086319675 44.62675158031099R的列向量为相应的近似特征向量R =0.804598963384540.570

14、432516844010.16500682363931-0.580702349614620.697759448610210.419424049176060.12411804571692-0.433288005375390.89266803187144（4）实验体会：本次程序是由老师提供的，但是同学们都认真阅读过了，我发现此程序是相当复杂的， 要对矩阵进行迭代等很多操作，我还了解到我们在处理较大的问题是必须要使用到程序，因 此让我对程序产生的浓厚的兴趣，同时我们也认识到我们所编写的程序稳定性很差，因此我 们还需要更多的练习。实验四插值法4.1实验目的掌握插值法的基本思路和步骤；培养编程与上机调试

15、能力。4.2实验要求用Matlab和插值中的某种具体算法编写代码并执行，完成解决具体问题。4.3实验内容计算书上的习题Matlab Spline Tools题目：x =00.20000.40000.60000.8000y =1.00001.22141.49181.82212.2255x0 =0.1500原理：设计思想：对应程序：function s=Newton(x,y,x0,nn)%nn为Newton插值多项式的次数nx=length(x);ny=length(y);if nx=nywarning (向量x与y的长度应该相同)return;endm=length(x0);for i=1:my

16、y=y;t=0;kk=1;while(kk=nn)kk=kk+1;for k=kk:nxyy(k) = (yy(k)-yy(kk-1)/(x(k)-x(kk-1);endendt=yy(1);for k=2:nn+1u=1.0;j=1；while(j Newton(x,y,x0,3)ans =1.1619 Newton(x,y,x0,4)ans =1.1618实验体会：本实验通过对牛顿插值公式的程序化，是我们知道了可以用插值公式来解决函数问题， 虽然我们只练习了牛顿公式，但是我们对插值公式有了很好的理解，对于在实验过程中遇到 的问题都在老师的细心辅导下得到了很好的解决，本次试验受益匪浅。实验五

17、最小二乘法5.1实验目的掌握最小二乘法的基本思路和拟合步骤；培养编程与上机调试能力。5.2实验要求用Matlab和插值中的某种具体算法编写代码并执行，完成解决具体问题。5.3实验内容计算书上的习题题目：已知一组实验数据如下:X1234Y1.953.053.553.85求问题的最小二乘法。原理：已知数据对q, y )( j = 1,2,n )，求多项式p(x)=祝 axi (m n)ii =0使得(a ,a，,a ) = (lLaxi - y 为最小，这就是一个最小二乘问题。0 i n = 11。i j jJ设计思想：用线性函数p(x) = a + bx为例，拟合给定数据(x , y.), i

18、= 1,2,m。算法描述：Z步骤1：输入m值，及(x., y)i = 1,2,m。步骤2：建立法方程组ATAX = AY。步骤3：解法方程组。步骤4：输出p(x) = a + bx。对应程序：x=1:1:4;y=1.95,3.05,3.55,3.85;p=Polyfit(x,y,3);x1=1:0.1:4;y1=polyval(p,x1);plot(x,y,*r,x1,y1,b)实验结果：ans =0.0667-0.70002.7333-0.1500图形(如果可视化)最小二乘拟合的专门的函数，我们可以通过调用相应的函数就可以达到我们需要的拟合结果，但要了解最小二乘拟合的思想，所以老师建议我们自

19、己编写相应的代码，通过本次试验， 也极大的激发了我们的编程兴趣，是我们认识到了我们的不足之处，我们编写的程序都是在 正常的运行情况下基本不出错，但是但出现异常时整个程序就会崩溃，也就是说我们编写的 程序不够健壮，这需要我们在以后的学习中不断进取。实验六复化求积公式与龙贝格算法6.1实验目的掌握复化求积公式与龙贝格算法的基本思路和迭代步骤;培养编程与上机调试能力。6.2实验要求编写程序，要求给出相应习题的实验结果。6.3实验内容题目：用龙贝格算法计算：I = j1业dx0 x（2）原理：龙贝格算法利用外推法，提高了计算精度，加快了收敛速度。R = R+ w1ki，ji ,k = 2,3,k, j

20、k, j-14 j-1 一 1对每一个k, j从2做到k，一直做到R -RRk 1 = T （f）（复化梯度求积公式），h =（3）k ,k rf,k-11小于给定的精度是停止计算。其中 k , 2k-1=T (f )hk设务思想：步骤1：输入区间端点a, b，精度控制值e，循环次数M，定义函数f (x)，取n = 1， h = b - a; R= h (f (a) + f (b)/2步骤 2： for k = 2 to M(Rk -1,1+h -1Vj = 2 to(Rk ,1forif瑚.=RR 一 Rk ,kk -1,k-1笑2 f (a +(2i-1)h )i=1kk,j-1k,j-1

21、 e/2-R )/ (4 j-1-1)k-1, j-1退出循环步骤3：数据积分近似值Rk ,k（4）对应程序：function s=romberg1(fun,a,b,ep)%龙贝格算法-求积公式%此算法只外推到Romberg序列if nargin=3ep=1.0e-6;elseif nargin=eparea=0.0;n=n+1;h=h/2;for i=1:marea=area+feval(fun,h*(2*i-1)+a);endt(n+1,1)=0.5*t(n,1)+area*h;m=2*m;if n4for j=1:3for i=1:n-jt(i,j+1) = (4(j)*t(i+1,j)

22、-t(i,j)/(4(j)-1);endendt1=t(n-4,4);t2=t(n-3,4);endenddisp(用Romberg序列求得积分近似值为)；s=t2;endfunction y=fun(x) if x=0y=1;y=sin(x)/xromberg1(fun,0,1,10e-6)实验结果：ans =0.94608307036726实验体会：以前不知道用程序怎么求积分，但通过本次实验我知道了利用程序是可以求积分的，它 主要用的的有限次的迭代，将连续的无限多点的问题有线化，然后就可以利用计算机来解决 积分问题，本程序由于比较复杂，所以老师给出了，虽然程序是由老师给出的，但是我们在 调

23、用该程序时也基本明白了程序实现的，在实验过程中我们仔细阅读了程序，所以基本上可 以再现，但我个人觉得这还远远不够，我们需要自己独立思考学习编写，不断的提高自己的 编程能力。实验七常微分方程数值解法7.1实验目的掌握常微分方程数值解的常用算法；培养编程与上机调试能力.求解常微分方程初值、边值问题的解.7.2实验题目：下述方面的相应习题（1）用改进的欧拉公式，求解常微分方程初值问题的解（2）用四阶龙格-库塔公式解初值问题（3）求解常微分方程边值问题的解7.3实验要求（1）选择一种计算机语言设计出改进欧拉法和四阶龙格-库塔法方法求解常微分方程初 值问题的程序，观察运行结果.（2）利用Matlab求解

24、常微分方程初值问题函数dsolve（）用于求解微分方程.Dy表示：dy/dt（t为缺省的自变量），Dny表示y对t 的n阶导数.Matlab6.5环境下操作如下： y=dsolve（Dy=y*y,y（0）=1）%求解题目 1 y=dsolve（Dy=y/t,y（2.0）=1）%求解题目 2（3）ode23， ode45（1）利用最小二乘法拟合通过改进欧拉法求出微分方程的一系列数值解的近似函数方程.并利用Matlab的绘图功能画出函数的曲线。（1）题目1：y =x2+y2,0 x1,y（0）=0,h=0.1题目2：y =3*y/（1+x）,0 x1,y（0）=1,h=0.2（2）原理：（3）设计

25、思想：（4）对应程序1：function s=Eulerc（fun,x0,xn,y0,h）%h为步长，改进的求微分方程数值解的Euler公式。if nargin5errorreturnendn=(xn-x0)/h;for i=1:nt1=y0+h*feval(fun,x0,y0);x0=x0+h;t2=y0+h*feval(fun,x0,t1);y0=(t1+t2)/2;y(i)=y0;ends=y;return对应程序2：function s=RK(fun,x0,xn,y0,h)%h为步长if nargin Eulerc(t2,x0,xn,y0,0.1)2.60003.5290ans =1.

26、11101.25151.43611.68802.04885.371510.348338.1343实验结果2: RK(t3,x0,xn,y0,0.2)ans =（6）图形（如果可视化）01.00000.20001.72750.40002.74300.60004.09420.80005.82921.00007.9960（7）实验体会通过本次试验，我们领略到了程序的魅力，是我们对程序产生了浓厚的兴趣，虽然在实 验过程中我们遇到了比想象中要大的问题，但最后通过老师的悉心辅导和我们自己的不懈努 力，最终都得到了圆满的解答，本次实验让我们看到了利用程序来求解微分方程的可能，我 感觉到到程序也是一种艺术。实

27、验八非线性方程求根8.1实验目的：掌握二分法、牛顿迭代法等常用的非线性方程迭代算法;培养编程与上机调试能力.8.2实验题目及参考结果：题目求方程/(x) x3 + x2 - 3x -3 = 0在1.5附近的根.原方程的根为x 1.7320518.3实验要求：(1)选择一种计算机语言设计出二分法和牛顿法的程序，并且选择不同的初值，观察所需 的迭代次数和迭代结果； (2)分析二分法和牛顿法在非线性方程求根中的优缺点和收敛速度二分法简单易行，但只有线性收敛，且仅限于求实根；牛顿法也是一种简单的迭代法，具有二阶收敛速度(在单根邻近处)的特点，但对初值 的选择比较苛刻，否则可能不收敛.Matlab6.1

28、环境下操作如下：f=1 1 -3 -3;roots(f)%多项式求根fx=x3+x2-3*x-3;fzero(fx,1.5)%非线性方程式的实根(1)实验题目：求方程/(x) - x3 + x2 3x 3 0在1.5附近的根.原方程的根为x 1.732051(2)基本思想及原理: 二分法 对所给的非线性方程进行初步确定其根所在区间为”,仞，此时有f (a)* f (b) 0是否成立，成立则将a x，否则令b x，对此利 02000用该方法，令x0逐步逼近方程的解x*。牛顿迭代法 作切线，则该切线与x轴的交点为，当I xk+ xk 10 disp（该区间无根！） return;endwhile abs(a-b)=epsfor k=1:mc=(a+b)/2if myf(a)*myf(c)=eps for k=1:m x0=x1;x1=x0-dao(x0) k=k+1endenddisp(x1);function y=dao(x) y=(x3) + (x2)-(3*x)-3; y1=3*(x”2)+2*x-3;y=y/y1;实验结果：二分法 midpart(1,2,0.00001,10) c =1.50000000000000