数值代数实验报告

本科实验报告课程名称：计算机数值方法实验项目：计算机数值方法实验实验地点：专业班级：学号：学生姓名：xxx指导教师：xxx太原理工大学学生实验报告学院名称软件学院专业班级1217班学号201200xxxx 学生姓名xx 实验日期2014.05.21 成绩课程名称数值计算方法实验题目实验一方程求解一、实验目的和要求熟悉使用、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。

选择上述方法中的两种方法求方程：二分法f(x)=x3+4x2-10=0在[1,2]内的一个实根，且要求满足精度|x*-x n|<0.5×10-5二、主要设备笔记本 HP ProBook 6470b 一台编译软件：VC++6.0三、实验内容和原理函数f(x)在区间（x，y）上连续，先在区间（x，y）确定a与b，若f(a)，f(b)异号，说明在区间(a，b)内存在零点，然后求f[(a+b)/2]。

假设F(a)<0,F(b)>0,a<b，①如果f[(a+b)/2]=0，该点即为零点；②如果f[(a+b)/2]<0，则区间（(a+b)/2，b)内存在零点，(a+b)/2≥a；③如果f[(a+b)/2]>0，则区间(a,(a+b)/2)内存在零点，(a+b)/2≤b；返回①重新循环，不断接近零点。

通过每次把f(x)的零点所在区间收缩一半的方法，使区间内的两个端点逐步逼近函数零点，最终求得零点近似值。

四、操作方法与实验步骤1. 二分法：#include<stdio.h>#include<stdlib.h>#include<math.h>int main(){double a=1.0, b=2.0;double x,s;printf(" An\t\tBn\t\tF(Xn)\n");while(1){x=(a+b)/2;s=pow(x,3)+4*x*x-10;if (-0.000005 < s && s < 0.000005){break;}else if(s < 0){a=x;}else if(s > 0){b=x;}printf("%f\t%f\t%f\n",a,b,s);}printf("X的值为：%f\n",x);printf("误差：\t%f\n",s);return 0;}2. 割线法：#include"stdio.h"#include"math.h"int main(){float c,a=1.0,b=2.0;printf("每次得到的X的近似值：\n");while(1){c=b-(b*b*b+4*b*b-10)*(b-a)/(b*b*b+4*b*b-(a*a*a+4*a*a));if(fabs(b-c)<0.5*0.00001)break;b=c;printf("%f\n",b);}printf("X的值为：%f\n",c);}五、实验结果与分析二分法割线法分析：由程序知，使用二分法和割线法均能计算出方程的根，但利用割线法要比二分法计算的次数少，并且能够较早的达到精度要求。

数值代数实验报告数值代数实验报告引言：数值代数是一门研究数值计算方法和算法的学科，它在科学计算和工程应用中起着重要的作用。

本实验报告旨在通过实际的数值计算问题，探讨数值代数的应用和效果。

实验一：线性方程组求解线性方程组求解是数值代数中的一个重要问题。

在实验中，我们使用了高斯消元法和LU分解法两种求解线性方程组的方法，并对比了它们的效果。

首先，我们考虑一个3×3的线性方程组：2x + 3y - z = 54x - 2y + 2z = 1x + y + z = 3通过高斯消元法，我们将该方程组转化为上三角形式，并得到解x=1, y=2, z=0。

而通过LU分解法，我们将该方程组分解为LU两个矩阵的乘积，并得到相同的解。

接下来，我们考虑一个更大的线性方程组，例如10×10的方程组。

通过比较高斯消元法和LU分解法的运行时间，我们可以发现LU分解法在处理大规模方程组时更加高效。

实验二：特征值与特征向量计算特征值与特征向量计算是数值代数中的另一个重要问题。

在实验中，我们使用了幂法和QR方法两种求解特征值与特征向量的方法，并对比了它们的效果。

首先，我们考虑一个3×3的矩阵：1 2 34 5 67 8 9通过幂法，我们可以得到该矩阵的最大特征值为15.372，对应的特征向量为[0.384, 0.707, 0.577]。

而通过QR方法，我们也可以得到相同的结果。

接下来，我们考虑一个更大的矩阵，例如10×10的矩阵。

通过比较幂法和QR 方法的运行时间，我们可以发现QR方法在处理大规模矩阵时更加高效。

实验三：奇异值分解奇异值分解是数值代数中的一种重要技术，它可以将一个矩阵分解为三个矩阵的乘积，从而实现数据降维和信息提取的目的。

在实验中，我们使用了奇异值分解方法，并通过实际的数据集进行了验证。

我们选取了一个包含1000个样本和20个特征的数据集，通过奇异值分解，我们将该数据集分解为三个矩阵U、S和V的乘积。

数值实验报告二一、实验名称解线性方程组的列主元素高斯消去法和LU 分解法二、实验目的通过数值实验，从中体会解线性方程组选主元的必要性和LU 分解法的优点，以及方程组系数矩阵和右端向量的微小变化对解向量的影响。

三、实验内容解下列两个线性方程组（1） ⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛--11134.981.4987.023.116.427.199.103.601.3321x x x （2） ⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----15900001.582012151526099999.23107104321x x x x 四、算法描述1、 列主元素高斯消去法记： ij ij a a =1)( (i, j = 1,2,3n )i i b b =1)( (i = 1,2,3n )消元过程：对于k = 1,2,3n(1) 选行号k i ，使)()(max k i ni k k k i k k a a ≤≤=。

(2) 交换)(k kj a 与)(k j i k a （j = k, k+1,k+2n ）以及)()(k i k k k b b 与所含的数值。

（3）对于i = k, k+1,k+2n ,计算)()(k kkk ik ik a a m =)()()1(k kj ik k ij k ij a m a a -=+ （j = k, k+1,k+2n ）)()()1(k k ik k i k i b m b b -=+回代过程：)(n nnn n a b x = )()1)()(/(k kk j n k j k kj k k k a x a a x ∑+=-= （k = n-1, n-2, n-3 1 ）在此算法中的)(k k i k a 称为第k 个列主元素，它的数值总要被交换到第k 个主对角线元素的位置上。

2、 LU 分解法通过MATLAB 自有的函数，把系数矩阵A 分解成A=LU ，其中：L 是下三角矩阵，U 是上三角矩阵，这时方程组Ax=b 就可以分解成两个容易求解的三角形方程组Ly=b ，Ux=y 。

第1篇一、实验目的1. 理解数值计算的基本概念和常用算法；2. 掌握Python编程语言进行数值计算的基本操作；3. 熟悉科学计算库NumPy和SciPy的使用；4. 分析算法的数值稳定性和误差分析。

二、实验内容1. 实验环境操作系统：Windows 10编程语言：Python 3.8科学计算库：NumPy 1.19.2，SciPy 1.5.02. 实验步骤（1）Python编程基础1）变量与数据类型2）运算符与表达式3）控制流4）函数与模块（2）NumPy库1）数组的创建与操作2）数组运算3）矩阵运算（3）SciPy库1）求解线性方程组2）插值与拟合3）数值积分（4）误差分析1）舍入误差2）截断误差3）数值稳定性三、实验结果与分析1. 实验一：Python编程基础（1）变量与数据类型通过实验，掌握了Python中变量与数据类型的定义方法，包括整数、浮点数、字符串、列表、元组、字典和集合等。

（2）运算符与表达式实验验证了Python中的算术运算、关系运算、逻辑运算等运算符，并学习了如何使用表达式进行计算。

（3）控制流实验学习了if-else、for、while等控制流语句，掌握了条件判断、循环控制等编程技巧。

（4）函数与模块实验介绍了Python中函数的定义、调用、参数传递和返回值，并学习了如何使用模块进行代码复用。

2. 实验二：NumPy库（1）数组的创建与操作通过实验，掌握了NumPy数组的基本操作，包括创建数组、索引、切片、排序等。

（2）数组运算实验验证了NumPy数组在数学运算方面的优势，包括加、减、乘、除、幂运算等。

（3）矩阵运算实验学习了NumPy中矩阵的创建、操作和运算，包括矩阵乘法、求逆、行列式等。

3. 实验三：SciPy库（1）求解线性方程组实验使用了SciPy库中的线性代数模块，通过高斯消元法、LU分解等方法求解线性方程组。

（2）插值与拟合实验使用了SciPy库中的插值和拟合模块，实现了对数据的插值和拟合，并分析了拟合效果。

数值计算方法实验报告实验目的：通过实验验证不同数值计算方法在求解数学问题时的精度和效率，并分析其优缺点。

实验原理：实验内容：本实验选取了三个典型的数值计算问题，并分别采用了二分法、牛顿迭代法和梯度下降法进行求解。

具体问题和求解方法如下：1. 问题一：求解方程sin(x)=0的解。

-二分法：利用函数值的符号变化将解空间不断缩小，直到找到满足精度要求的解。

-牛顿迭代法：通过使用函数的斜率来逼近方程的解，并不断逼近真实解。

-梯度下降法：将方程转化为一个极小化问题，并利用梯度下降的方式逼近极小值点，进而找到方程的解。

2.问题二：求解函数f(x)=x^2-3x+2的极小值点。

-二分法：通过确定函数在一个区间内的变化趋势，将极小值所在的区间不断缩小，从而找到极小值点。

-牛顿迭代法：通过使用函数的导数和二阶导数来逼近极小值点，并不断逼近真实解。

-梯度下降法：将函数转化为一个极小化问题，并利用梯度下降的方式逼近极小值点，进而找到函数的极小值点。

3. 问题三：求解微分方程dy/dx = -0.1*y的解。

-二分法：通过离散化微分方程，将微分方程转化为一个差分方程，然后通过迭代计算不同点的函数值，从而得到函数的近似解。

-牛顿迭代法：将微分方程转化为一个积分方程，并通过迭代计算得到不同点的函数值，从而得到函数的近似解。

-梯度下降法：将微分方程转化为一个极小化问题，并利用梯度下降的方式逼近极小值点，从而得到函数的近似解。

实验步骤：1.编写代码实现各个数值计算方法的求解过程。

2.对每个数值计算问题，设置合适的初始值和终止条件。

3.运行程序，记录求解过程中的迭代次数和每次迭代的结果。

4.比较不同数值计算方法的精度和效率，并分析其优缺点。

实验结果：经过实验测试，得到了如下结果：-问题一的二分法迭代次数为10次，求解结果为x=0；牛顿迭代法迭代次数为4次，求解结果为x=0；梯度下降法迭代次数为6次，求解结果为x=0。

-问题二的二分法迭代次数为10次，求解结果为x=1;牛顿迭代法迭代次数为3次，求解结果为x=1;梯度下降法迭代次数为4次，求解结果为x=1-问题三的二分法迭代次数为100次，求解结果为y=e^(-0.1x);牛顿迭代法迭代次数为5次，求解结果为y=e^(-0.1x);梯度下降法迭代次数为10次，求解结果为y=e^(-0.1x)。

试验项目名称：平方根法与改进平方根法实验内容：先用你熟悉的计算机语言将平方根法和改进平方根法编写成通用的子程序，然后用你编写的程序求解对称正定方程组Ax=b，其中，A=[101 10 1…1 10 11 10]100*100b随机生成，比较计算结果，评论方法优劣。

实验要求：平方根法与改进的平方根的解法步骤；存储单元，变量名称说明；系数矩阵与右端项的生成；结果分析。

实验报告姓名：罗胜利班级：信息与计算科学0802 学号：u200810087实验一、平方根法与改进平方根法先用你所熟悉的计算机语言将平方根法和改进的平方根法编成通用的子程序，然后用你编写的程序求解对称正定方程组AX=b，其中系数矩阵为40阶Hilbert矩阵，即系数矩阵A的第i行第j列元素为=，向量b的第i个分量为=.平方根法函数程序如下：function [x,b]=pingfanggenfa(A,b)n=size(A);n=n(1);x=A^-1*b; %矩阵求解disp('Matlab自带解即为x');for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endend %Cholesky分解for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n); %前代法A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1); %回代法disp('平方根法的解即为b');endfunction [x]=ave(A,b,n) %用改进平方根法求解Ax=b L=zeros(n,n); %L为n*n矩阵D=diag(n,0); %D为n*n的主对角矩阵S=L*D;for i=1:n %L的主对角元素均为1L(i,i)=1;for i=1:nfor j=1:n %验证A是否为对称正定矩阵if (eig(A)<=0)|(A(i,j)~=A(j,i)) %A的特征值小于0或A非对称时，输出wrong disp('wrong');break;endendendD(1,1)=A(1,1); %将A分解使得A=LDL Tfor i=2:nfor j=1:i-1S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);endD(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endy=zeros(n,1); % x，y为n*1阶矩阵x=zeros(n,1);for i=1:ny(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i); %通过LDy=b解得y的值endfor i=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n)); %通过L T x=y解得x的值改进平方根法函数程序如下：function b=gaijinpinfanggenfa(A,b)n=size(A);n=n(1);v=zeros(n,1);for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);end %LDL'分解B=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;EndA=tril(A); %得到L和Dfor j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n); %前代法A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1); %回代法disp('改进平方根法解得的解即为b');end调用函数解题：clear;clc;n=input('请输入矩阵维数：');b=zeros(n,1);A=zeros(n);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);b(i)=b(i)+1/(i+j-1);endend %生成hilbert矩阵[x,b]=pingfanggenfa(A,b)b=gaijinpinfanggenfa(A,b)运行结果：请输入矩阵维数：40Warning: Matrix is close to singular or badly scaled.Results may be inaccurate. RCOND = 6.570692e-020.> In pingfanggenfa at 4In qiujie at 10Matlab自带解即为x平方根法的解即为bx =1.60358.96850.85621.01950.9375-50.2500-3.0000-16.000024.0000-49.5000-30.000039.000022.0000-64.0000-12.00002.000010.2500-10.5000-1.0000-10.875083.000046.0000-98.000012.0000-69.000068.000021.0000-50.7188-8.7500-8.0000 112.00006.0000 -68.750022.000044.0000 -28.00008.0000 -44.000012.0000b =1.0e+007 *0.0000-0.00000.0001-0.0004-0.00140.0424-0.29801.1419-2.73354.2539-4.30182.7733-1.19890.5406-0.36880.3285-0.44380.4621-0.25130.05650.0000-0.00510.0071-0.0027-0.0031-0.00190.00090.0002-0.0002-0.00060.00040.0001-0.00020.00010.0000-0.00000.0000-0.0000-0.0000改进平方根法解得的解即为bb =1.0e+024 *0.0000-0.00000.0001-0.00120.0139-0.09540.4208-1.21012.0624-1.0394-3.33436.2567-0.2463-7.45942.80303.69900.7277-1.7484-0.4854-3.60100.25325.18621.44100.8738-4.56541.04224.0920-2.7764-2.2148-0.89530.36654.89671.04160.1281-4.3387-1.1902-2.83348.4610-3.6008实验二、利用QR分解解线性方程组：利用QR分解解线性方程组Ax=b，其中A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12];b=[32 26 38 30];求解程序如下：定义house函数：function [v,B]=house(x)n=length(x);y=norm(x,inf);x=x/y;Q=x(2:n)'*x(2:n);v(1)=1;v(2:n)=x(2:n);if n==1B=0;elsea=sqrt(x(1)^2+Q);if x(1)<=0v(1)=x(1)-a;elsev(1)=-Q/(x(1)+a);endB=2*v(1)^2/(Q+v(1)^2);v=v/v(1);endend进行QR分解：clear;clc;A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12];b=[32 26 38 30];b=b';x=size(A);m=x(1);n=x(2);d=zeros(n,1);for j=1:n[v,B]=house(A(j:m,j));A(j:m,j:n)=(eye(m-j+1)-B*(v')*v)*A(j:m,j:n);d(j)=B;if j<mA(j+1:m,j)=v(2:m-j+1);endend %QR分解R=triu(A); %得到R D=A;I=eye(m,n);Q=I;for i=1:nD(i,i)=1;endH=tril(D);M=H';for i=1:nN=I-d(i)*H(1:m,i)*M(i,1:m);Q=Q*N;end %得到Qb=(Q')*b; %Q是正交阵for j=n:-1:2b(j)=b(j)/R(j,j);b(1:j-1)=b(1:j-1)-b(j)*R(1:j-1,j);endb(1)=b(1)/R(1,1); %回带法运行结果如下：R =18.7617 9.8072 15.7769 11.08640 9.9909 9.3358 7.53410 0 5.9945 9.80130 0 0 -0.5126Q =0.8528 -0.4368 -0.2297 -0.17090.2132 0.7916 -0.4594 -0.34170.4264 0.3822 0.2844 0.76890.2132 0.1911 0.8095 -0.5126b=1.000000000000001.000000000000010.9999999999999881.00000000000001实验三、Newton下山法解非线性方程组：3x-cos（yz）-=0，-81+sinz+1.06=0，exp（-xy）+20z+=0；要求满足数值解=满足或.定义所求方程组的函数：Newtonfun.mfunction F = Newtonfun(X)F(1,1)=3*X(1)-cos(X(2)*X(3))-1/2;F(2,1)=X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;F(3,1)=exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3;End向量求导：Xiangliangqiudao.mfunction J=xiangliangqiudao()syms x y zX=[x,y,z];F=[3*X(1)-cos(X(2)*X(3))-1/2;X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;exp(-X(1)*X(2))+20*X(3)+(10 *pi-3)/3];J=jacobian(F,[x y z]);End代值函数：Jacobi.mfunction F=Jacobi(x)F=[ 3,x(3)*sin(x(2)*x(3)), x(2)*sin(x(2)*x(3));2*x(1), -162*x(2)-81/5,cos(x(3));-x(2)/exp(x(1)*x(2)),-x(1)/exp(x(1)*x(2)),20];End方程组求解：format long; %数据表示为双精度型X1=[0,0,0]';eps=10^(-8);k=1;i=1;X2=X1-Jacobi(X1)^(-1)*Newtonfun(X1);while (norm(subs(X2-X1,pi,3.1415926),2)>=eps)&&(norm(Newtonfun(X1),2)>=eps) if norm(Newtonfun(X2),2)<norm(Newtonfun(X1),2) %判断先后两次迭代的大小X1=X2;B=inv(Jacobi(X2));C=Newtonfun(X2);X2=X2-B*C;i=i+1;elsev=1/(2^k); %引入下山因子X1=X2;B=inv(Jacobi(X2));C=Newtonfun(X2);X2=X2-v*B*C;k=k+1;endendj=i+k-1 %迭代次数X=X2 %输出结果运行结果如下：j =5X =0.500000000000000 -0.000000000000000 -0.523598775598299。

一、引言很多情况下在工程中抽象出来的数学方程组是超定的，没有精确解，这样就需要找一个在某种意义下最接近精确解的解。

设A 是m*n 的实矩阵，22A (A )T x bQ x b -=-,这样2min A y b -就等价与2min (A )T Q x b -,为方便求解，需要A T Q 是上三角矩阵，这样引入QR 分解就比较求解方便。

二、豪斯霍尔德变换豪斯霍尔德变换（Householder transformation ）又称初等反射（Elementary reflection ），最初由A.C Aitken 在1932年提出，Alston Scott Householder 在1958年指出了这一变换在数值线性代数上的意义。

这一变换将一个向量变换为由一个超平面反射的镜像，是一种线性变换。

其变换矩阵被称作豪斯霍尔德矩阵，在一般内积空间中的类比被称作豪斯霍尔德算子，超平面的法向量被称作豪斯霍尔德向量。

三、理论依据householder 变换：任意的一个向量x ，经过一个正交变换2**T H I W W =- 后总可以变为一个与之范数相等的另一个向量Hx 。

如上图中所示，记v Hx x =-，2/w v v =，2**TH I W W =-上述H 既为所要求的householder 变换。

具体操作时将需要变化的矩阵的每一列当做一个向量，第一列变为除了第一个元素都为0的向量…第i 列变为除了前i 个元素都为0的向量...为了保证在每次变化时不改变已经变好的0元素，第i 次只变化每列第i 个元素到第n 个元素，每个变化矩阵的形式是[I 0;0 H]。

算法如下：qr_Houserhoilder.m function [Q,R]=qr_Houserhoilder(A) [m,n]=size(A); Q=eye(n); for i=1:nu=house(A(i:m,i));P=eye(m-i+1)-2/(norm(u)^2)*u*u'; A(i:m,i:n)=P*A(i:m,i:n);PP=blkdiag(eye(i-1),P);Q=Q*PP;endR=triu(A);程序house.mfunction u=house(Ai)Ai(1)=Ai(1)+sign(Ai(1))*norm(Ai);u=Ai/norm(Ai);四、数值实验结果：为了保证在一定阶数下，所做的数值试验具有代表性，我们随机产生3个n阶矩阵，随着n的不断变化，得到如下一系列图：（1）A=rand(4,4)A =0.3027 0.7894 0.1793 0.39220.2364 0.7134 0.4350 0.39360.4690 0.0742 0.9747 0.11150.5381 0.6768 0.7245 0.8650[Q,R]=qr_Houserhoilder(A)Q =-0.3735 -0.5369 0.4602 0.6003-0.2916 -0.5439 -0.7840 -0.0668-0.5786 0.6445 -0.2679 0.4219-0.6639 -0.0208 0.3190 -0.6761R =-0.8106 -0.9951 -1.2387 -0.90010 -0.7781 0.2803 -0.37070 0 -0.2886 0.11800 0 0 -0.3286test_qr3err =1.0e-10 *0.1950 0.0120 0.0119 0.3029-111 1.52 2.53 3.54（2）A=rand(5,5)A =0.7575 0.7005 0.6420 0.8233 0.88770.6373 0.4698 0.5815 0.5290 0.37700.6643 0.3036 0.9458 0.7210 0.12950.1056 0.0254 0.3903 0.1048 0.30270.2813 0.3787 0.5302 0.0627 0.1156[Q,R]=qr_Houserhoilder(A)Q =-0.6161 0.4378 0.2968 -0.5805 -0.0598-0.5183 -0.0381 0.2932 0.7098 -0.3743-0.5403 -0.6729 -0.2440 -0.1030 0.4303-0.0859 -0.1847 -0.5513 -0.2507 -0.7693-0.2288 0.5657 -0.6801 0.2928 0.2817R =-1.2295 -0.9280 -1.3629 -1.1944 -0.86480 0.2940 -0.1496 -0.1288 0.29660 0 -0.4455 0.1231 0.09690 0 0 -0.1847 -0.30310 0 0 0 -0.3388test_qr3err =1.0e-10 *0.1722 0.0119 0.0120 0.4239-111 1.52 2.53 3.54（3）A=rand(6,6)A =0.5961 0.8422 0.4279 0.3473 0.2731 0.34860.1335 0.2860 0.2157 0.6181 0.6306 0.19300.4479 0.0467 0.3284 0.4588 0.4228 0.95480.9430 0.7296 0.7951 0.0920 0.0322 0.31140.7581 0.7811 0.4338 0.1512 0.6185 0.12710.8647 0.1359 0.1735 0.9036 0.4393 0.8964 [Q,R]=qr_Houserhoilder(A)Q =-0.3571 -0.5571 0.2402 0.2038 0.4576 0.5035-0.0800 -0.2530 -0.1880 0.8253 -0.1645 -0.4314-0.2684 0.3532 -0.4622 0.2419 -0.3487 0.6399-0.5650 -0.0874 -0.6012 -0.3372 0.2898 -0.3376-0.4542 -0.3266 0.3151 -0.2560 -0.7185 -0.0773-0.5181 0.6217 0.4823 0.1991 0.2045 -0.1760R =-1.6690 -1.1737 -0.9944 -0.8854 -0.7882 -1.09430 -0.7594 -0.2803 0.3166 -0.0940 0.58290 0 -0.3472 0.1832 0.1390 -0.10880 0 0 0.8020 0.5966 0.50220 0 0 0 -0.4714 -0.02290 0 0 0 0 0.4305test_qr3err =1.0e-10 *0.1518 0.0119 0.0120 0.2845-111 1.52 2.53 3.54五、总结与分析由上面的数值试验结果和绘图可知：比较之下，我自己写的houserhoilder程序误差是最大的，CGS与MGS的误差相差不大，系统自带的matlab软件速度是最快的，我自己写的houserhoilder程序还是最慢的。

数值代数实验报告Numerical Linear Algebra And ItsApplications学生所在学院：理学院学生所在班级：计算数学10-1学生姓名：戈东潮指导教师：于春肖教务处2012年12月实验一实验名称： Poisson 方程边值问题的五点差分格式 实验时间: 2012年12月13日 星期四 实验成绩： 一、实验目的：通过上机利用Matlab 数学软件实现Poisson 方程的边值问题的五点差分格式的线性代数方程组来认真解读Poisson 方程边值问题的具体思想与方法，使我们掌握得更加深刻，并且做到学习与使用并存，增加学习的实际动手性，不再让学习局限于书本和纸上，而是利用计算机学习来增加我们的学习兴趣。

二、实验内容：利用Poisson 方程来解下列问题：⎪⎩⎪⎨⎧∂∈=∈=∂∂+∂∂-ΩΩ),(,),(),(,sin sin )(y x y x u y x y x y u x u 0222222πππ 其中}{的边界是ΩΩ∂<<∈Ω,1,0),(y x y x 。

边值问题的解释是y x y x u ππsin sin ),(=（1）用例题中模型问题的方法取N=5,10（也可以取N=20）列出五点差分格式的线性代数方程组三、实验过程：系数矩阵A 的Matlab 实现函数如下：%文件名：poisson.m function T=poisson(N) for i=1:N a(i)=4; end b=diag(a); for j=1:N-1b(j,j+1)=-1;endfor j=2:Nb(j,j-1)=-1;endfor i=1:NI(i)=-1;endI=diag(I);for i=1:N:N*Nfor j=1:NT(i+j-1,i:i+N-1)=b(j,:);endendfor i=1:N:N*N-Nfor j=1:NT(i+j-1,i+N:i+N-1+N)=I(j,:);endendfor i=1+N:N:N*Nfor j=1:NT(i+j-1,i-N:i-1)=I(j,:);endend求解常数项b的Matlab实现函数如下：%函数名：constant.mfunction b=constant(N)n=N+1;h=1/n;i=1;for y=1/n:1/n:N/nfor x=1/n:1/n:N/nf(i)=2*pi*pi*sin(pi*x)*sin(pi*y);i=i+1;endendb=h*h*f;b=b';四、实验结果(总结/方案)在Matlab运行窗口输入>>A=poisson(5) Enter输出结果A =Columns 1 through 114 -1 0 0 0 -1 0 0 0 0 0-1 4 -1 0 0 0 -1 0 0 0 00 -1 4 -1 0 0 0 -1 0 0 00 0 -1 4 -1 0 0 0 -1 0 00 0 0 -1 4 0 0 0 0 -1 0-1 0 0 0 0 4 -1 0 0 0 -10 -1 0 0 0 -1 4 -1 0 0 00 0 -1 0 0 0 -1 4 -1 0 00 0 0 -1 0 0 0 -1 4 -1 00 0 0 0 -1 0 0 0 -1 4 00 0 0 0 0 -1 0 0 0 0 40 0 0 0 0 0 -1 0 0 0 -10 0 0 0 0 0 0 -1 0 0 00 0 0 0 0 0 0 0 -1 0 00 0 0 0 0 0 0 0 0 -1 00 0 0 0 0 0 0 0 0 0 -10 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 Columns 12 through 220 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 00 -1 0 0 0 0 0 0 0 0 00 0 -1 0 0 0 0 0 0 0 00 0 0 -1 0 0 0 0 0 0 0 -1 0 0 0 -1 0 0 0 0 0 04 -1 0 0 0 -1 0 0 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 00 -1 4 -1 0 0 0 -1 0 0 00 0 -1 4 0 0 0 0 -1 0 00 0 0 0 4 -1 0 0 0 -1 0 -1 0 0 0 -1 4 -1 0 0 0 -10 -1 0 0 0 -1 4 -1 0 0 00 0 -1 0 0 0 -1 4 -1 0 00 0 0 -1 0 0 0 -1 4 0 00 0 0 0 -1 0 0 0 0 4 -10 0 0 0 0 -1 0 0 0 -1 40 0 0 0 0 0 -1 0 0 0 -10 0 0 0 0 0 0 -1 0 0 00 0 0 0 0 0 0 0 -1 0 0 Columns 23 through 250 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 0 4 -1 0 -1 4 -10 -1 4 在运行窗口输入 >>b=constant(5) Enter 常数项的输出结果b =[ 0.1371 0.2374 0.2742 0.2374 0.1371 0.2374 0.4112 0.4749 0.4112 0.2374 0.2742 0.4749 0.5483 0.4749 0.2742 0.2374 0.4112 0.4749 0.4112 0.2374 0.1371 0.2374 0.2742 0.2374 0.1371]T所以A*u=b 的五点差分格式为j i j i j i j i j i j i f h u u u u u ,,,,,,211114=-----+-+ （i,j=1,2……）实验二实验名称： 用Jacobi 迭代法和SOR 迭代法求解方程组 实验时间: 2012年12月13日 星期四 实验成绩： 一、实验目的：通过上机利用Matlab 数学软件采用Jacobi 迭代法和SOR 迭代法来求解方程组，求解过程中了解两种迭代方法的基本思想与迭代过程，并分析两种迭代方法的收敛性和实用用以及他们的异同点。

数值分析实验报告模板篇一：数值分析实验报告(一)(完整)数值分析实验报告12345篇二：数值分析实验报告实验报告一题目：非线性方程求解摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。

本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。

利用二分法求解给定非线性方程的根，在给定的范围内，假设f(x,y)在[a,b]上连续，f(a)xf(b) 直接影响迭代的次数甚至迭代的收敛与发散。

即若x0 偏离所求根较远，Newton法可能发散的结论。

并且本实验中还利用利用改进的Newton法求解同样的方程，且将结果与Newton法的结果比较分析。

前言：（目的和意义）掌握二分法与Newton法的基本原理和应用。

掌握二分法的原理，验证二分法，在选对有根区间的前提下，必是收敛，但精度不够。

熟悉Matlab语言编程，学习编程要点。

体会Newton使用时的优点，和局部收敛性，而在初值选取不当时，会发散。

数学原理：对于一个非线性方程的数值解法很多。

在此介绍两种最常见的方法：二分法和Newton法。

对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b) Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式xk?1?xk?f(xk) f'(xk)产生逼近解x*的迭代数列{xk}，这就是Newton法的思想。

当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。

另外，若将该迭代公式改进为xk?1?xk?rf(xk) 'f(xk)其中r为要求的方程的根的重数，这就是改进的Newton 法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。

程序设计：本实验采用Matlab的M文件编写。

其中待求解的方程写成function的方式，如下function y=f(x);y=-x*x-sin(x);写成如上形式即可，下面给出主程序。

数值代数实验数值线性代数实验一一、实验名称：矩阵的LU分解.二、实验目的：用不选主元的LU分解和列主元LU分解求解线性方程组Ax=b, 并比较这两种方法.三、实验内容与要求（1）用所熟悉的计算机语言将不选主元和列主元LU分解编成通用的子程序，然后用编写的程序求解下面的84阶方程组将计算结果与方程组的精确解进行比较，并就此谈谈你对Gauss消去法的看法.（2）写出追赶法求解三对角方程组的过程，并编写程序求该实验中的方程组Gauss消去法：用消去法解方程组的基本思想是用逐次消去未知数的方法把原来方程组Ax=b化为与其等价的三角方程组，而求解三角方程组就容易了。

换句话说，上述过程就是用行的初等变换将原方程组系数矩阵化为简单形式，从而将求解原方程组的问题转化为求解简单方程组的问题。

利用Gauss消去法对线性方程组Ax=b进行求解。

用MATLAB建立m文件DelGauss.m,程序如下：function x=DelGauss(a,b)[n,m]=size(a);nb=length(b);det=1;x=zeros(n,1);for k=1:n-1for i=k+1:nif a(k,k)==0returnendm=a(i,k)/a(k,k);for j=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);for k=n:-1:1for j=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);End在matlab中输入如下：结果如下：方程组的精确解为x1=x2=…=x84=1.0000，与Gauss消去法求得的解差距很大，所得结果不够准确，计算简单但其消元过程有时不能进行到底而使求解出现解失真的情况。

数值线性代数实验二一、实验名称：实对称正定矩阵的A的Cholesky分解.二、实验目的：用平方根法和改进的平方根方法求解线性方程组Ax=b.三、实验内容与要求用所熟悉的计算机语言将Cholesky分解和改进的Cholesky分解编成通用的子程序，然后用编写的程序求解对称正定方程组Ax=b，其中（1）b随机的选取，系数矩阵为100阶矩阵（2）系数矩阵为40阶Hilbert矩阵，即系数矩阵A的第i行第j列元素为，向量b的第i个分量为（3）用实验一的程序求解这两个方程组，并比较所有的计算结果，然后评价各个方法的优劣。

课程实验报告COURSE PAPER课程名称：数值代数与计算方法课程号：08305114授课教师：学号：姓名：所属：计算机科学与工程打印时间：2015评语：题目一：算法：1、对于问题一：2、对于问题二：直接编写递归函数程序，算出三个差分方程的10个近似值程序:1）%main.mclc;clear all;a=1;b=-2;c=-3;[x1,x2]=roots(a,b,c)%roots.mfunction [x1,x2]=roots(a,b,c)d=sqrt(b*b-4*a*c);if d>0x1=(-2*c)/(b+d);x2=(-b-d)/(2*a);elseif d<0x1=(-b+d)/(2*a);x2=(-2*c)/(b-d);end2）%solu.mfunction [X,R,P,Q]=solu(X0,R0,P0,P1,Q0,Q1) X(1)=X0;R(1)=R0;P(1)=P0;P(2)=P1;Q(1)=Q0;Q(2)=Q1;for i=1:9X(i+1)=X(i)/2;X(i)=X(i+1);endfor i=1:9R(i+1)=R(i)/2;R(i)=R(i+1);endfor i=3:10P(i)=3/2*P(i-1)-1/2*P(i-2);P(i-1)=P(i);P(i-2)=P(i-1);endfor i=3:10Q(i)=5/2*Q(i-1)-Q(i-2);Q(i-1)=Q(i);Q(i-2)=Q(i-1);end%Xn.mclc;clear all;X0=1;R0=0.994;P0=1;P1=0.497;Q0=1;Q1=0.497;[X,R,P,Q]=solu(X0,R0,P0,P1,Q0,Q1);for i=1:10x(i)=i;endplot(x,X-R,'r*');hold on结果：1）x1 =3x2 =-12）题目二：1、此题我们以第一个式子为例，求出g(Xn+1)= (2/(2+3*Xn-Xn.^3))^(1/2)2、利用二分法：3、利用牛顿迭代法，并根据题目提示的立方根算法：此题我们以第三个式子为例：，设A=71/3 ,求出x的3次方便是近似值。

数值线性代数实验题目：数值线性代数专业：信息与计算科学班级：班姓名：山东科技大学2013年 1 月16日实验报告说明学院：信息学院专业：信息班级10-2 姓名：一、主要参考资料：（1）《Matlab数值计算-案例分析》北京航空出版（2）《Matlab数值分析》机械工业出版二、课程设计应解决的主要问题：（1）平方根（2）QR方法（3）最小二乘法三、应用软件：（1）Matlab7.0（2）数学公式编辑器四、发出日期：课程设计完成日期：指导教师签字：系主任签字：指导教师对课程设计的评语指导教师签字：年月日一、问题描述先用你所熟悉的计算机语言将平方根和改进的平方根法编成写通用的子程序，然后用你编写的程序求解对称正定方程组b x =A ，其中 （1）b 随机的选取，系数矩阵位100阶矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡1011101110111011101110（2）系数矩阵为40阶Hilbert 矩阵，即系数矩阵A 的第i 行第j 列元素为11-+=j i a ij ，向量b 的第i 个分量为∑=-+=nj i j i b 111。

二、分析与程序1. 平方根法函数程序如下：function [x,b]=pingfanggenfa(A,b) n=size(A); n=n(1);x=A^-1*b; disp('Matlab 自带解即为x'); for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k); for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k); endend for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('平方根法的解即为b');endfunction [x]=ave(A,b,n)求解Ax=bL=zeros(n,n);D=diag(n,0);S=L*D;for i=1:n %L的主对角元素均为1L(i,i)=1;endfor i=1:nfor j=1:nif (eig(A)<=0)|(A(i,j)~=A(j,i))disp('wrong');break;endendendD(1,1)=A(1,1);for i=2:nfor j=1:i-1S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);endD(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endy=zeros(n,1);x=zeros(n,1);for i=1:ny(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i); endfor i=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n));end2.改进平方根法函数程序如下：function b=gaijinpinfanggenfa(A,b)n=size(A);n=n(1);v=zeros(n,1);for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);end %LDL'分解B=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;EndA=tril(A);for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j); endb(n)=b(n)/A(n,n);A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('改进平方根法解得的解即为b'); end3.调用函数解题：clear;clc;n=input('请输入矩阵维数：');b=zeros(n,1); A=zeros(n);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);b(i)=b(i)+1/(i+j-1);endend[x,b]=pingfanggenfa(A,b)b=gaijinpinfanggenfa(A,b)4.运行结果：请输入矩阵维数：40Warning: Matrix is close to singular or badly scaled.Results may be inaccurate. RCOND = 6.570692e-020. > In pingfanggenfa at 4In qiujie at 10Matlab自带解即为x平方根法的解即为bx =1.60358.96850.85621.01950.9375-50.2500-3.0000-16.000024.0000-49.5000-30.000039.000022.0000-64.0000 -12.00002.000010.2500 -10.5000-1.0000 -10.875083.000046.0000 -98.000012.0000 -69.000068.000021.000017.0000 -50.7188-8.7500-8.0000 112.00006.0000 -68.750022.000044.0000 -28.00008.0000 -44.000012.0000b =1.0e+007 *0.0000-0.00000.0001-0.0004-0.00140.0424-0.29801.1419-2.73354.2539-4.30182.7733-1.19890.5406-0.36880.3285-0.44380.4621-0.25130.05650.0000-0.00510.0071-0.0027-0.00310.0036-0.00190.00090.0002-0.0002-0.00060.00040.0001-0.00020.00010.0000-0.00000.0000-0.0000-0.0000改进平方根法解得的解即为bb =1.0e+024 *0.0000-0.00000.0001-0.00120.0139-0.09540.4208-1.21012.0624-1.0394-3.33436.2567-0.2463-7.45942.80303.69900.7277-1.7484-0.4854-3.60100.25325.1862-2.12991.44100.8738-4.56541.04224.0920-2.7764-2.2148-0.89530.36654.89671.04160.1281 -4.3387 -1.1902 -2.8334 8.4610 -3.6008一、问题描述先用你所熟悉的计算机语言将算法2.5.1编成写通用的子程序，然后用你编写的程序完成下面两个计算任务：（1） 估计5到20阶Hilbert 矩阵的∞范数条件数；（2）设n *n n R 11-1-1-111-1-101-1001A ∈⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=先随机的选取n R x ∈，并计算出x A b n =；然后再用列主元Gauss 消去法求解该方程组，假定计算解为xˆ。

暨南大学本科实验报告专用纸课程名称数值代数成绩评定实验项目名称(填写实验所属章节名称) 指导教师刘娟实验项目编号0701******* 实验项目类型验证实验地点机房学生姓名王伟文学号2010051727学院信息科学技术学院数学系信息与计算科学专业2010级实验时间2011年9月 1 日～12月30日温度24℃第一章1.实验选题：（写出你在该项目所选的实验题目）第一章上机练习12.谈谈你对该算法的理解：（简单谈一下你是如何理解该算法的？）对算法的理解：先将84阶的矩阵A分解为一个下三角矩阵L和上三角矩阵U,先考虑下三角形方程组Ly=b,利用前代法解出y,再考虑上三角形方程组Ux=y;利用回代法解出x。

列主元gauss消去法Ax b PA LU Ly Pb Ux y=⇔===,,;3.实验内容（将实验程序及其实验结果粘贴，最好对程序各部分注释清楚，比如设置了哪些函数，这些函数的输入输出是什么，具有什么功能？）实验程序对矩阵A进行LU分解的程序function [ L,U ] = LUfac( A )for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endL=tril(A,0);for i=1:nL(i,i)=1;endU=triu(A,0);end利用前代法解出y值的程序function [ b ] = TSL( L,b )n=size(L,1);for j=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);end利用回代法解出x值的程序function [ b ] = TSU( U,b )n=size(U,1);for j=n:-1:2b(j)=b(j)/U(j,j);b(1:j-1)=b(1:j-1)-b(j)*U(1:j-1,j);endb(1)=b(1)/U(1,1);end主函数程序（‘生成84阶的矩阵A’）A=eye(84);A=6*A;for i=2:84A(i,i-1)=8;A(i-1,i)=1;end（‘生成84乘1的矩阵b’）b=ones(84,1);b=b*15;b(1)=7;b(84)=14;[L,U]=LUfac(A);（‘调用函数LUfac对矩阵A进行分解’）y=TSL(L,b);（‘调用函数TSL求解Ly=b方程’）x=TSU(U,y);（‘调用函数TSU求解UX=b方程’）用MATLAB求解得到的结果x’ans =1.0e+008 *Columns 1 through 70.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 8 through 140.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 15 through 210.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 22 through 280.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 29 through 350.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 36 through 420.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 43 through 490.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 50 through 560.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 Columns 57 through 63-0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 Columns 64 through 700.0000 -0.0000 0.0000 -0.0001 0.0002 -0.0003 0.0007 Columns 71 through 77-0.0013 0.0026 -0.0052 0.0105 -0.0209 0.0419 -0.0836 Columns 78 through 840.1665 -0.3303 0.6501 -1.2582 2.3487 -4.0263 5.3684列主元gauss消去法function [L,U,P]=Lufac(A)n=size(A,1);P=eye(n,n);for i=1:n-1[r,m]=max(abs(A(i:n,i)));m=m+i-1;A([i,m],:)=A([m,i],:);P([i,m])=P([m,i]);if A(i,i)~=0A(i+1:n,i)=A(i+1:n,i)/ A(i,i);A(i+1:n,i+1:n)= A(i+1:n,i+1:n)-A(i+1:n,i)*A(i,i+1:n); endendU=triu(A);L=tril(A,-1)+eye(n,n);endA=eye(84);A=6*A;for i=2:84A(i,i-1)=8; A(i-1,i)=1; endb=ones(84,1);b=b*15;b(1)=7;b(84)=14;[L,U,P]=Lufac(A); b=P*b;y=TSL(L,b);x=TSU(U,y);求解结果为x =1.0e+025 *0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00010.0002-0.00040.0007-0.00150.0030-0.00600.0119-0.02380.0475-0.09460.1878-0.36960.7154-1.33552.2894-3.05251.教师评语、评分：（请认真对待每次数值代数实验，期末交实验报告，将计算实验成绩）2.实验选题：第三章上机习题1（3）。

数值代数上机报告一、 目的意义： 把矩阵A 分解成一个下三角阵与一个上三角阵的乘积，即 A=LR ，其中L 为下三角阵，R 为上三角阵，这样原线性方程组就可以化为⎩⎨⎧==y Rx b Ly 的求解问题，方便求解。

二、 算法：1) 输入系数矩阵 A ；2) 利用公式iji k ij ij ij r l a r ∑-=-=11和ii ki i k jk ji ji r r l a l /)(11∑-=-=交错进行，计算得出矩阵L 和R ；3) 回带到⎩⎨⎧==y Rx b Ly 中得出原线性方程组的解。

三、 源程序：#include <stdio.h>#include <string.h>#include <math.h>#include <stdlib.h>#define N 100main(){int i,j,k,s,n;printf("请输入系数矩阵A 的阶数：n= ");scanf("%d",&n);floata[N][N]={0},L[N][N]={0},R[N][N]={0},sigma1,sigma2,b[N],y[N],x[N];/*为L 主对角线元素赋1*/for(i=0;i<n;i++){L[i][i]=1;}printf("请输入系数矩阵A ：\n"); /*输入系数矩阵A*/报告1 Doolittle 分解for(i=0;i<n;i++){for(j=0;j<n;j++)scanf("%f",&a[i][j]);}for(k=0;k<n;k++){for(j=k;j<n;j++) /*计算矩阵R*/{sigma1=0;for(s=0;s<=k-1;s++)sigma1+=L[k][s]*R[s][j];R[k][j]=a[k][j]-sigma1;}for(i=k;i<n;i++) /*计算矩阵L*/{sigma2=0;for(s=0;s<=k-1;s++)sigma2+=L[i][s]*R[s][k];L[i][k]=(a[i][k]-sigma2)/R[k][k];}}printf("\n A矩阵为：\n");/*输出矩阵L、R*/ for(i=0;i<n;i++){for(j=0;j<n;j++)printf("%5.1f ",a[i][j]);printf("\n");}printf("\n L矩阵为：\n");for(i=0;i<n;i++){for(j=0;j<n;j++)printf("%5.1f ",L[i][j]);printf("\n");}printf("\n R矩阵为：\n");for(i=0;i<n;i++){for(j=0;j<n;j++)printf("%5.1f ",R[i][j]);printf("\n");}printf("请输入b矩阵\n",i+1);for(i=0;i<n;i++)scanf("%f",&b[i]);for(i=0;i<n;i++)/*回代法求解方程组Ly=b*/{sigma1=0;for(k=0;k<=i-1;k++)sigma1+=L[i][k]*y[k];y[i]=b[i]-sigma1;}for(i=n-1;i>=0;i--){sigma2=0;for(k=i+1;k<n;k++)sigma2+=R[i][k]*x[k];x[i]=(y[i]-sigma2)/R[i][i];}printf("解得x为：\n");for(i=0;i<n;i++)printf("%5.1f ",x[i]);printf("\n");}四、计算结果与分析：分析：运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想五、 参考文献：[1]刑志栋. 矩阵数值分析. 陕西： 陕西科学技术出版社， 2005[2]谭浩强. C 语言程序设计. 北京：清华大学出版社，2005一、 目的意义： 对称正定矩阵是实践中经常遇到的一种特殊矩阵类型矩阵，由于矩阵本身的流量好兴致，使得cholesky 分解在存储和运算量上较一般消去法节省一半左右，且解的精度高，cholesky 分解方法是目前计算机求解该类问题最有效的方法之一。

数值代数课程设计实验报告姓名： 班级： 学号： 实验日期：一、实验名称 代数的数值解法 二、实验环境实验一、平方根法与改良平方根法一、实验要求：用熟悉的运算机语言将不选主元和列主元Gasuss 消元法编写成通用的子程序，然后用编写的程序求解以下方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--⨯14151515157681681681681681612321n n n n n x x x x x x 用所编的程序别离求解40、84、120阶方程组的解。

二、算法描述及实验步骤GAuss 程序如下：（1）求A 的三角分解：LU A =;（2）求解b y =L 得y ； （3）求解y x =U 得x ；列主元Gasuss 消元法程序如下： 1求A 的列主元分解：LU PA =;2求解b y P L =得y ；3求解y x U 得x ；三、调试进程及实验结果：%----------------方程系数---------------->> A1=Sanduijiaozhen(8,6,1,40); >> A2=Sanduijiaozhen(8,6,1,84); >> A3=Sanduijiaozhen(8,6,1,120); >> b1(1)=7;b2(1)=7;b3(1)=7; >> for i=2:39 b1(i)=15; end>> b1(40)=14; >> for i=2:83 b2(i)=15; end>> b2(40)=14; >> for i=2:119 b1(i)=15; end>> b3(120)=14;%----------------方程解---------------->> x11=GAuss(A1,b1') >> x12=GAuss Zhu (A1,b1') >> x21=GAuss(A2,b2') >> x22=GAuss Zhu (A3,b3') >> x31=GAuss(A3,b3') >> x32=GAuss Zhu (A3,b3') 运行结果：（n=40） GAuss 消元法的解即为x11 =列主元GAuss消元法的解即为x12 =1 1 1 1 1 1 1 1 1 1 111111111111111111111111111111六、源程序：function A=Sanduijiaozhen(a,b,c,n)%生成n阶以a,b,c为元素的三对角阵A=diag(b*ones(1,n),0)+diag(c*ones(1,n-1),1)+diag(a*ones(1,n-1),-1);function x=GAuss(A,b)n=length(b);x=b;%-------分解---------------for i=1:n-1for j=i+1:nmi=A(j,i)/A(i,i);b(j)=b(j)-mi*b(i);for k=i:nA(j,k)=A(j,k)-mi*A(i,k);endAB=[A,b];endend%-----------回代------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1s=0;for j=i+1:ns=s+A(i,j)*x(j);endx(i)=(b(i)-s)/A(i,i);endfunction x=GAussZhu(A,b)n=length(b);x=b;%----------------------选主元---------------------for k=1:n-1a_max=0;for i=k:nif abs(A(i,k))>a_maxa_max=abs(A(i,k));r=i;endendif r>kfor j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%--------------消元-----------------for i=k+1:nm=A(i,k)/A(k,k);for j=k:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);endendif abs(A(n,n))==0return;endAbZhu=[A,b];%----------------回代-----------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1for j=i+1:nb(i)=b(i)-A(i,j)*x(j);endx(i)=b(i)/A(i,i);end实验二、平方根法与改良平方根法一、实验要求：用运算机语言将平方根法和改良的平方根法编成通用的子程序，然后用编写的程序求解对称正定方程组100阶方程组AX=b，二、算法描述及实验步骤：平方根法函数程序如下：一、求A 的Cholesky 分解：L L A T=; 二、求解b y =L 得y ；3、求解y x =TL 得x ；改良平方根法函数程序如下：一、求A 的Cholesky 分解：T=LDL A ; 二、求解b y =L 得y ； 3、求解y x =TDL 得x ；三、调试进程及实验结果：clear;clc;%----------------方程系数---------------->> A=Sanduijiaozhen(1,10,1,100); >> b(1)=11; >> for i=2:99 b(i)=12; end>> b(100)=11;>> x1=Cholesky(A,b') >> x2=GJCholesky(A,b')运行结果：平方根法的解即为 x1 =改良平方根法解得的解即为x2 =四、源程序：function x=Cholesky(A,b)n=size(A);n=n(1);% x=A^-1*b;% disp('Matlab自带解即为x');%-----------------Cholesky分解-------------------for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endend%------------------前代法求解Ly=b----------------------------for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);%-----------------回代法求解L'x=y-----------------------------A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('平方根法的解即为');function b=GJCholesky(A,b)n=size(A);n=n(1);v=zeros(n,1);%----------------------LDL'分解-----------------------------for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);endB=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;end%-------------------前代法---------------------------A=tril(A); %取得L和Dfor j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);%-----------------回代法-----------------------------A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('改良平方根法解得的解即为');实验三、二次多项式拟合一、实验要求：用运算机语言编制利用QR分解求解线性最小二乘问题的通用子程序，用编写的程序求解一个二次多项式使在残向量的范数最小的意义下拟合下面的数据t-1 0iyi二、算法描述及实验步骤：QR分解求解程序如下：一、求A 的QR 分解; 二、计算b c 11T=Q ;3、求解上三角方程1c x =R 得x ；三、调试进程及实验结果：>> t=[-1 0 ]; >> y=[ ]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图');运行后屏幕显示数据的散点图（略）.（3）编写以下MATLAB 程序计算)(x f 在),(i i y x 处的函数值，即输入程序 >> syms a b c >> t=[-1 0 ]; >> fi=a.*t.^2+ b.*t+c%运行后屏幕显示关于 ,,a b c 的线性方程组fi =[a-b+c,9/16*a-3/4*b+c,1/4*a-1/2*b+c,c,1/16*a+1/4*b+c,1/4*a+1/2*b+c,9/16*a+3/4*b +c]编写构造残向量2范数的MATLAB 程序>> y=[ ]; >> y=[ ];>> fy=fi-y; fy2=fy.^2; J=sum(fy.^2); 运行后屏幕显示误差平方和如下 J=(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2为求,,a b c 使J 达到最小，只需利用极值的必要条件0J a ∂=∂，0J b ∂=∂，0Jc∂=∂，取得关于,,a b c 的线性方程组，这能够由下面的MATLAB 程序完成，即输入程序 >> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3) 运行后屏幕显示J 别离对,,a b c 的偏导数如下 Ja11 =451/128*a-63/32*b+43/8*c-887/128 Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8解线性方程组112131000Ja Ja Ja ===，，，输入以下程序 >> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8]; >> C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f 及其系数C 如下 C =f =924/2999*x^2+10301/11996*x+4204/2999 故所求的拟合曲线为2()0.30810.8581 1.4018f x x x =++四、源程序：>> t=[-1 0 ]; >> y=[ ]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图'); >> syms a b c>> t=[-1 0 ]; >> fi=a.*t.^2+ b.*t+c fi =[ a-b+c, 9/16*a-3/4*b+c, 1/4*a-1/2*b+c, c, 1/16*a+1/4*b+c, 1/4*a+1/2*b+c, 9/16*a+3/4*b+c]>> y=[ ]; >> y=[ ];>> fy=fi-y; fy2=fy.^2; J=sum(fy.^2) J =(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2>> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3)Ja11 =451/128*a-63/32*b+43/8*c-887/128Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8>> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8];>> C=B/A, f=poly2sym(C)C =f =924/2999*x^2+10301/11996*x+4204/2999>>。

一、 实验目的与要求计算曲边梯形面积dx e I x⎰=20，可用内接梯形的面积n S 来进行逼近，用动画形象地描述这一逼近的过程（如果],[)(b a C x f ∈，按照定积分的定义，有∑⎰-=+=≈=10)()(n i n b a h ih a f S dx x f I ，其中n a b h -=。

这实际上是一个最简单的复化求积公式）实验要求：（1）用数值实验说明I S n n =∞→lim （其中2,,)(===b o a e x f x ）；（2）设计一副动画，模拟I S n →的过程一、 实验方案(程序源文件)a=input('请输入积分上限a= ');b=input('请输入积分上限b= ');n=[10,50,100,200,10000];m=length(n);for k=1:ms=0;i=1;h=(b-a)/n(k);g=exp(a);while i<n(k)y=exp(a+(i-1)*h);g=[g;y];s=s+y*h;i=i+1;endx=0:0.1:2;plot(x,exp(x),'r')hold onf=sym('2.71828^x');fprintf('曲边梯形面积')digits(5)I=int(f,'x',a,b)title(['y=exp(x),n=',num2str(n(k))])fprintf('内接梯形的面积')sx2=0.5*h:h:2-0.5*h;bar(x2,g,1.0,'c')pause(2)hold offend二、实验结果和数据处理(运行结果) 结果：请输入积分上限a= 0请输入积分上限b= 2曲边梯形面积I =6.3891内接梯形的面积s = 4.5615曲边梯形面积I = 6.3891内接梯形的面积s = 5.9782曲边梯形面积I =6.3891内接梯形的面积s = 6.1805曲边梯形面积I =6.3891内接梯形的面积s = 6.2840曲边梯形面积I =6.3891内接梯形的面积s = 6.3869曲边梯形面积I =6.3891内接梯形的面积s = 4.5615曲边梯形面积I =6.3891内接梯形的面积s = 5.9782曲边梯形面积I =6.3891内接梯形的面积s = 6.1805曲边梯形面积I =6.3891内接梯形的面积s = 6.2840曲边梯形面积I =6.3891内接梯形的面积s = 6.3869二、结论依图可知,当n趋于无穷大时，内接阶梯形面积趋近于曲边梯形的面积。

数值线性代数第二版徐树方高立张平文上机习题第一章实验报告（供参考）上机习题1.先用你所熟悉的的计算机语言将不选主元和列主元Gauss 消去法编写成通用的子程序；然后用你编写的程序求解84阶方程组；最后将你的计算结果与方程的精确解进行比较，并就此谈谈你对Gauss 消去法的看法。

Sol ：（1）先用matlab 将不选主元和列主元Gauss 消去法编写成通用的子程序，得到P U L ,,：不选主元Gauss 消去法：[])(,A GaussLA U L =得到U L ,满足LU A =列主元Gauss 消去法：[])(,,A GaussCol P U L =得到P U L ,,满足LU PA =（2）用前代法解()Pb or b Ly =，得y用回代法解y Ux =，得x求解程序为()P U L b A Gauss x ,,,,=（P 可缺省，缺省时默认为单位矩阵）（3）计算脚本为ex1_1代码%算法function [L,U]=GaussLA(A)n=length(A);for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endU=triu(A);L=tril(A);L=L-diag(diag(L))+diag(ones(1,n));end%算法function [L,U,P]=GaussCol(A)n=length(A);for k=1:n-1[s,t]=max(abs(A(k:n,k)));p=t+k-1;temp=A(k,1:n);A(k,1:n)=A(p,1:n);A(p,1:n)=temp;u(k)=p;if A(k,k)~=0A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); else break;endendL=tril(A);U=triu(A);L=L-diag(diag(L))+diag(ones(1,n));P=eye(n);for i=1:n-1temp=P(i,:);P(i,:)=P(u(i),:);P(u(i),:)=temp;end。