1、 山东建筑大学数值计算 A实验报告二级学院： 理学院 专 业： 信息与计算科学 指导教师： XXX 班级学号： 信计 XX 20081212XX 姓 名： XXXX 实验一 高斯消去法（1）题目：用选主元素法和高斯消去法求解下列方程组：已知方程增广矩阵：a =2 -1 0 0 6-1 -3 -2 0 1-1 3 -2 0 00 0 -3 5 1（2）原理： (1)高斯消去法：相对于约当消去法而言，总的来说就是将增广矩阵化为下三角阵。(2) 顺序选主元素法:相对于高斯消去法的唯一不同是要先按当前要排列元素所在列大小进行排列。（3）设计思想： (1)高斯消去法：首先让每一行的元素除以该行的主对角线
2、元素。然后利用此行使位于下一行主对角线以前的元素变为 0，依次类推。(2)选主元素法：在高斯消去法的基础上，每次进行化上三角阵之前，重新排列各方程的位置。（4）对应程序：（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(
3、主元素为零，消去法无法继续) ;break;else a(i,j)=a(i,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=2 ep=0.000001endm,n=size(a);if m=ndisp(输入错误，系数矩证阵只能是方阵)endif n=length(b)disp(输入错
4、误,常数项的个数应与方程的个数相同)endfor k=1:n-1p=a(k,k);I=k;for i=k:nif abs(a(i,k)abs(p)p=a(i,k);I=i;endendif pa =2 -1 0 0 -1 -3 -2 0 -1 3 -2 0 0 0 -3 5 b =6 1 0 1 gauss1(a,b)方程组的解为ans =35/12 -1/6 -41/24 -33/40 （6）实验体会：主元消去法和高斯消去法的确是两个非常锻炼人编程的方法，在编写程序时，需要使用的大量的循环和分支结构，但无论是高斯消去法还是高斯列主元法，它们的原理还算不难理解，通过变成能够较好的理解它们。实验
5、二 解线性方程组的迭代法2.1 实验目的 掌握解线性方程组的雅可比迭代和高斯-塞德尔迭代算法； 培养编程与上机调试能力.2.2 实验要求：（1）选择一种计算机语言（Matlab）设计出雅可比（Jacobi）Gauss-Seidel、SOR 迭代法，迭代法的算法程序，并且选择不同的迭代次数，观察输出结果；（2）利用 Matlab 求方程组2.3 实验内容计算书上的习题（1）题目： 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)%Jac
6、obi 迭代法求解方程组if nargin=3ep=1.0e-6;elseif nargin=ep)I=i;endendif p=eps)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)2-a(p,q)*sin(2*theta);G(p,q)=(a(q,q)-a(p,p)*(cos(theta)*(sin(th
7、eta)+a(p,q)*(cos(2*theta);G(q,p)=G(p,q);for 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)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+gasi2);G
8、(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)*(1/sqrt(1
9、+t2);if (i=p)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.87999781719956j =3lanmda =43.87999781719956u1 =0.185867901919500.4460
10、32888095101.00000000000000经典 Jacobi 法：G 的对角线即为近似的特征值G =0.03838874652912 0.00000000000020 0.000000000137200.00000000000020 2.33485967315988 -0.000000863196750.00000000013720 -0.00000086319675 44.62675158031099R 的列向量为相应的近似特征向量R =0.80459896338454 0.57043251684401 0.16500682363931-0.58070234961462 0.697
11、75944861021 0.419424049176060.12411804571692 -0.43328800537539 0.89266803187144（4）实验体会：本次程序是由老师提供的，但是同学们都认真阅读过了，我发现此程序是相当复杂的，要对矩阵进行迭代等很多操作，我还了解到我们在处理较大的问题是必须要使用到程序，因此让我对程序产生的浓厚的兴趣，同时我们也认识到我们所编写的程序稳定性很差，因此我们还需要更多的练习。实验四 插值法4.1 实验目的 掌握插值法的基本思路和步骤； 培养编程与上机调试能力。4.2 实验要求用 Matlab 和插值中的某种具体算法编写代码并执行，完成解决具体
12、问题。4.3 实验内容计算书上的习题Matlab Spline Tools（1）题目： x = 0 0.2000 0.4000 0.6000 0.8000y = 1.0000 1.2214 1.4918 1.8221 2.2255x0 = 0.1500（2）原理：（3）设计思想： （4）对应程序：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:myy=y;t=0;k
13、k=1;while(kk Newton(x,y,x0,3)ans =1.1619 Newton(x,y,x0,4)ans =1.1618（7）实验体会：本实验通过对牛顿插值公式的程序化，是我们知道了可以用插值公式来解决函数问题，虽然我们只练习了牛顿公式，但是我们对插值公式有了很好的理解，对于在实验过程中遇到的问题都在老师的细心辅导下得到了很好的解决，本次试验受益匪浅。实验五 最小二乘法5.1 实验目的 掌握最小二乘法的基本思路和拟合步骤； 培养编程与上机调试能力。5.2 实验要求用 Matlab 和插值中的某种具体算法编写代码并执行，完成解决具体问题。5.3 实验内容计算书上的习题（1） 题目
14、： 已知一组实验数据如下：X 1 2 3 4Y 1.95 3.05 3.55 3.85求问题的最小二乘法。（2）原理：已知数据对 ，求多项式,1,2jxyn 0()()miipxan使得 为最小，这就是一个最小二乘问题。20110(,)nijjjiay（3）设计思想： 用线性函数 为例，拟合给定数据 。()pxb,1,2ixym算法描述：步骤 1：输入 值，及 。m,1,2iym步骤 2：建立法方程组 。TAXY步骤 3：解法方程组。步骤 4：输出 。()pxab（4）对应程序：x=1:1:4;y=1.95,3.05,3.55,3.85;p=Polyfit(x,y,3);x1=1:0.1:4;
15、y1=polyval(p,x1);plot(x,y,*r,x1,y1,b)（5）实验结果： ans =0.0667 -0.7000 2.7333 -0.1500（6）图形（如果可视化）（7）实验体会：通过本次试验，我们发现 MATLAB 工具箱的功能是很强大的，MATLAB 工具箱提供了最小二乘拟合的专门的函数，我们可以通过调用相应的函数就可以达到我们需要的拟合结果，但要了解最小二乘拟合的思想，所以老师建议我们自己编写相应的代码，通过本次试验，也极大的激发了我们的编程兴趣，是我们认识到了我们的不足之处，我们编写的程序都是在正常的运行情况下基本不出错，但是但出现异常时整个程序就会崩溃，也就是说我
16、们编写的程序不够健壮，这需要我们在以后的学习中不断进取。实验六 复化求积公式与龙贝格算法6.1 实验目的 掌握复化求积公式与龙贝格算法的基本思路和迭代步骤； 培养编程与上机调试能力。6.2 实验要求编写程序，要求给出相应习题的实验结果。6.3 实验内容（1）题目： 用龙贝格算法计算： 10sinxId（2）原理：龙贝格算法利用外推法，提高了计算精度，加快了收敛速度。 ,1,1, 2,3.4kjkjkjjR对每一个 从 2 做到 ，一直做到 小于给定的精度是停止计算。其中,kj ,1,kk（复化梯度求积公式） ，,1khRTf2bah（3）设计思想： 步骤 1：输入区间端点 ，精度控制值 ，循环
17、次数 ，定义函数 ，取,aeM()fx， n;hb1,()/Rf步骤 2：for to kM2,1,11 1,1,1, /2/4k ki jkjjkjkjfihfortRRife 退 出 循 环步骤 3：数据积分近似值 。,k（4）对应程序：function s=romberg1(fun,a,b,ep)% 龙贝格算法 -求积公式% 此算法只外推到 Romberg 序列if nargin=3ep=1.0e-6;elseif nargin=ep area=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.
18、5*t(n,1)+area*h;m=2*m;if n4 for j=1:3for i=1:n-jt(i,j+1)=(4(j)*t(i+1,j)-t(i,j)/(4(j)-1);endendt1=t(n-4,4);t2=t(n-3,4); endend disp(用 Romberg 序列求得积分近似值为);s=t2;endfunction y=fun(x)if x=0y=1;y=sin(x)/xromberg1(fun,0,1,10e-6)（5）实验结果： ans =0.94608307036726（6）实验体会：以前不知道用程序怎么求积分，但通过本次实验我知道了利用程序是可以求积分的，它主要用
19、的的有限次的迭代，将连续的无限多点的问题有线化，然后就可以利用计算机来解决积分问题，本程序由于比较复杂，所以老师给出了，虽然程序是由老师给出的，但是我们在调用该程序时也基本明白了程序实现的，在实验过程中我们仔细阅读了程序，所以基本上可以再现，但我个人觉得这还远远不够，我们需要自己独立思考学习编写，不断的提高自己的编程能力。实验七常微分方程数值解法7.1 实验目的 掌握常微分方程数值解的常用算法； 培养编程与上机调试能力.求解常微分方程初值、边值问题的解.7.2 实验题目：下述方面的相应习题（1）用改进的欧拉公式，求解常微分方程初值问题的解 （2）用四阶龙格-库塔公式解初值问题（3）求解常微分方
20、程边值问题的解7.3 实验要求（1） 选择一种计算机语言设计出改进欧拉法和四阶龙格-库塔法方法求解常微分方程初值问题的程序，观察运行结果.（2） 利用 Matlab 求解常微分方程初值问题函数 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） 利用最小二乘法拟合通过改进欧拉法求出微分方程的一系列数值解的近似函数方程.并利
21、用 Matlab 的绘图功能画出函数的曲线。（1）题目 1： y=x2+y2,0 Eulerc(t2,x0,xn,y0,0.1)ans =1.1110 1.2515 1.4361 1.6880 2.0488 2.6000 3.5290 5.3715 10.3483 38.1343实验结果 2: RK(t3,x0,xn,y0,0.2)ans =0 1.00000.2000 1.72750.4000 2.74300.6000 4.09420.8000 5.82921.0000 7.9960（6）图形（如果可视化）（7）实验体会通过本次试验，我们领略到了程序的魅力，是我们对程序产生了浓厚的兴趣，虽然
22、在实验过程中我们遇到了比想象中要大的问题，但最后通过老师的悉心辅导和我们自己的不懈努力，最终都得到了圆满的解答，本次实验让我们看到了利用程序来求解微分方程的可能，我感觉到到程序也是一种艺术。实验八 非线性方程求根8.1 实验目的： 掌握二分法、牛顿迭代法等常用的非线性方程迭代算法； 培养编程与上机调试能力.8.2 实验题目及参考结果：题目 求方程 在 1.5 附近的根.32()0fxx原方程的根为 1.758.3 实验要求：（1） 选择一种计算机语言设计出二分法和牛顿法的程序，并且选择不同的初值，观察所需的迭代次数和迭代结果;（2） 分析二分法和牛顿法在非线性方程求根中的优缺点和收敛速度二分法
23、简单易行，但只有线性收敛，且仅限于求实根；牛顿法也是一种简单的迭代法，具有二阶收敛速度（在单根邻近处）的特点，但对初值的选择比较苛刻，否则可能不收敛.(3) Matlab6.1 环境下操作如下：f=1 1 -3 -3;roots(f) %多项式求根fx=x3+x2-3*x-3;fzero(fx,1.5) %非线性方程式的实根（1）实验题目： 求方程 在 1.5 附近的根.32()0fxx原方程的根为 1.75（2）基本思想及原理：二分法对所给的非线性方程进行初步确定其根所在区间为 ，此时有 ，令,ba0)(*bfa，判断 是否成立，成立则将 ，否则令 ，对此利0bax0)(*xfa0xx用该方
24、法，令 逐步逼近方程的解 。0x*牛顿迭代法选取适当的初值 ，计算 ，然后过点 作切线，则该切线与 x 轴的交点)(0f )(,0fx为 ，据此得到迭代公式为 ，当 时，得到方)(00xf 1kk |1kx程的近似解。（3）对应程序：二分法function c=midpart(a,b,eps,m) format longformat compactif myf(a)*myf(b)0disp(该区间无根！)return;endwhile abs(a-b)=eps for k=1:m c=(a+b)/2 if myf(a)*myf(c)=epsfor k=1:mx0=x1;x1=x0-dao(x0
25、)k=k+1endenddisp(x1);function y=dao(x)y=(x3)+(x2)-(3*x)-3;y1=3*(x2)+2*x-3;y=y/y1;（4）实验结果： 二分法 midpart(1,2,0.00001,10)c =1.50000000000000k =2c =1.75000000000000k =3c =1.62500000000000k =4c =1.68750000000000k =5c =1.71875000000000k =6c =1.73437500000000k =7c =1.72656250000000k =8c =1.73046875000000k =
26、9c =1.73242187500000k =10c =1.73144531250000k =11c =1.73193359375000k =2c =1.73217773437500k =3c =1.73205566406250k =4c =1.73199462890625k =5c =1.73202514648438k =6c =1.73204040527344k =7c =1.73204803466797k =8c =1.73205184936523k =9c =1.73204994201660k =10c =1.73205089569092k =11输出解1.73205041885376
27、ans =1.73205041885376牛顿迭代法 newton(1.5,0.00001,10)x1 =1.73336066694000k =2x1 =1.73205192940947k =3x1 =1.73205080756970k =4x1 =1.73205080756888k =5x1 =1.73205080756888k =6x1 =1.73205080756888k =7x1 =1.73205080756888k =8x1 =1.73205080756888k =9x1 =1.73205080756888k =10x1 =1.73205080756888k =111.73205080756888ans =1.73205080756888（5）实验体会以前我们用手工来求解非线性方程的解时，我们都觉得比较繁琐，在编写程序时，虽然也比较繁琐，但是它有很好的通用性，正是以此编写到处可用，通过本次实验，我们感受到了程序的强大用途，使我们更渴望自己能够编写出好的程序，对于在实验过程中的问题，老师给我们精心的辅导，最终都得到了很好的解决，我们也从中写到了很多的知识。
