矩阵数值分析上机报告
- 格式:doc
- 大小:237.50 KB
- 文档页数:19
数值分析上机实践报告目录第一章标准题部分 (3)1.1 第一题 (3)1.1.1解题的理论依据、算法进行分析及应用条件 (3)1.1.2 MATLAB计算程序 (4)1..1.3 MATLAB运行结果 (7)1..1.4 结果分析与讨论 (7)1.2 第二题 (7)1.2.1解题的理论依据、算法进行分析及应用条件 (7)1.2.2 MATLAB计算程序 (8)1.2.2 MATLAB运行结果 (9)1.2.3 结果分析与讨论 (9)1.3 第三题 (9)1.3.1解题的理论依据、算法进行分析及应用条件 (9)1.3.2 MATLAB计算程序 (10)1.3.3 MATLAB运行结果 (11)1.3.4 结果分析与讨论 (11)1.4 第四题 (11)1.4.1解题的理论依据、算法进行分析及应用条件 (11)1.4.2 MATLAB计算程序 (12)1.4.3 MATLAB运行结果 (13)1.4.4 结果分析与讨论 (13)第二章自主题部分 (14)钢筋混凝土偏心受压柱配筋设计 (14)2.1 引言 (14)2.2 问题研究 (15)2.3 案例分析 (16)2.4 结论 (18)第一章 标准题部分1.1 第一题用三次样条插值的三弯矩法求f(-0.02)和f(2.56)。
1.1.1 解题的理论依据、算法进行分析及应用条件为使问题一般化,本题的三次样条求解使用了对于任意分化的三弯矩插值法。
若记各节点间距n x x x h j j j ,,2,111⋅⋅⋅=-=--,根据理论分析,三次样条插值函数为:))(6()()(6(6)(6)()(112112111131131-------------+--+-+-=j j j j j j j j j j j j jj j j h x x h M y h x x h M y h x x M h x x M x S其中),(1j j x x x -∈。
而式中的诸j M 可由如下的矩阵方程来确定:⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡N N d d d M M M10101111122212μλμλμ这里⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎨⎧-⋅⋅⋅=+=-=+=--'==---+='--=------+---+-)1,,2,1(1)](1[6),,(6)(6)(61111111111110000100N j h h h h h h y y h y h d x x x f h y y h y y h h d y h y y h d j j jj j j j j j N N N NN N j j j j j j j j j j j μλμ对于上述矩阵方程,由于其系数矩阵为一三对角阵,故我们用追赶法(Thomas 算法)来求解。
《矩阵分析与数值分析》实验报告院系:姓名:学号:所在班号:任课老师:一.设错误!未找到引用源。
,分别编制从小到大和从大到小的顺序程序计算错误!未找到引用源。
并指出有效位数。
程序如下:function sum3j=input('请输入求和个数 "j":');A=0;B=0;double B;double A;for n=2:jm=n^2-1;t=1./m;A=A+t;enddisp('从小到大:')s=Afor n=j:-1:2m=n^2-1;t=1./m;B=B+t;enddisp('从大到小:')s=B运行结果:>> sum3请输入求和个数 "j":100从小到大:s =0.740049504950495从大到小:s =0.740049504950495>> sum3请输入求和个数 "j":10000从小到大:s =0.749900004999506从大到小:s =0.749900004999500>> sum3请输入求和个数 "j":1000000从小到大:s =0.749999000000522从大到小:s =0.749999000000500二、解线性方程组1.分别Jacobi 迭代法和Gauss-Seidel 迭代法求解线性方程组。
⎪⎪⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----000121001210012100124321x x x x 迭代法计算停止的条件为:6)()1(3110max -+≤≤<-k j k j j x x 。
解:(1)Jacobi 迭代法程序代码: function jacobi(A, b, N) clc;clear;A=[-2 1 0 0;1 -2 1 0;0 1 -2 1;0 0 1 -2]; b=[-1 0 0 0]'; N=100;n = size(A,1); D = diag(diag(A)); L = tril(-A,-1); U = triu(-A,1); Tj = inv(D)*(L+U); cj = inv(D)*b; tol = 1e-06; k = 1;format longx = zeros(n,1); while k <= Nx(:,k+1) = Tj*x(:,k) + cj;disp(k); disp('x = ');disp(x(:,k+1)); if norm(x(:,k+1)-x(:,k)) < toldisp('The procedure was successful')disp('Condition ||x^(k+1) - x^(k)|| < tol was met after k iterations') disp(k); disp('x = ');disp(x(:,k+1)); break endk = k+1; end结果输出The procedure was successfulCondition ||x^(k+1) - x^(k)|| < tol was met after k iterations 60 x =0.799998799067310.599998427958700.399998056850090.19999902842505(2)Gauss-Seidel迭代法程序代码:function gauss_seidel(A, b, N)clc;clear;A=[-2 1 0 0;1 -2 1 0;0 1 -2 1;0 0 1 -2];b=[-1 0 0 0]';N=100;n = size(A,1);D = diag(diag(A));L = tril(-A,-1);U = triu(-A,1);Tg = inv(D-L)*U;cg = inv(D-L)*b;tol = 1e-06;k = 1;x = zeros(n,1);while k <= Nx(:,k+1) = Tg*x(:,k) + cg;disp(k); disp('x = ');disp(x(:,k+1));if norm(x(:,k+1)-x(:,k)) < toldisp('The procedure was successful')disp('Condition ||x^(k+1) - x^(k)|| < tol was met after k iterations') disp(k); disp('x = ');disp(x(:,k+1));breakendk = k+1;end结果输出The procedure was successfulCondition ||x^(k+1) - x^(k)|| < tol was met after k iterations31x =0.799999213979350.599998971085610.399999167590770.199999583795392. 用Gauss列主元消去法、QR方法求解如下方程组:⎪⎪⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛---017435222331325212214321x x x x (1)Gauss 列主元消去法 程序代码:function x=Gaussmain(A,b) clc;clear; format longA=[1 2 1 2;2 5 3 -2;-2 -2 3 5;1 3 2 3]; b=[4 7 -1 0]'; N=length(A); x=zeros(N,1); y=zeros(N,1); c=0; d=0;A(:,N+1)=b; for k=1:N-1 for i=k:4if c<abs(A(i,k))d=i;c=abs(A(i,k)); end endy=A(k,:);A(k,:)=A(d,:); A(d,:)=y; for i=k+1:N c=A(i,k);for j=1:N+1A(i,j)=A(i,j)-A(k,j)*c/A(k,k); end end endb=A(:,N+1);x(N)=b(N)/A(N,N); for k=N-1:-1:1x(k)=(b(k)-A(k,k+1:N)*x(k+1:N))/A(k,k); end结果输出 ans =18.00000000000000 -9.571428571428576.00000000000000-0.42857142857143(2)QR方法:程序代码function QR(A,b)clc;clear;format longA=[1 2 1 2;2 5 3 -2;-2 -2 3 5;1 3 2 3];b=[4 7 -1 0]';[Q,R]=qr(A)y=Q'*b;x=R\y结果输出Q =-0.31622776601684 -0.04116934847963 -0.75164602800283 0.57735026918962 -0.63245553203368 -0.49403218175557 -0.15032920560056 -0.57735026918963 0.63245553203368 -0.74104827263336 -0.22549380840085 -0.00000000000000 -0.31622776601684 -0.45286283327594 0.60131682240226 0.57735026918963 R =-3.16227766016838 -6.00832755431992 -0.94868329805051 2.84604989415154 0 -2.42899156029822 -4.65213637819829 -4.15810419644272 0 0 -0.67648142520255 -0.52615221960200 0 0 0 4.04145188432738 x =17.99999999999989-9.571428571428515.99999999999997-0.42857142857143三、非线性方程的迭代解法1.用Newton迭代法求方程()06cos22x=-++=-xexf x的根,计算停止的条件为:6110-+<-kkxx;编程如下:function newton(f,df,x,a,a0)syms xf=input('please enter your equation:') a0=input('please enter you x(0):');df=diff(f)e=1e-6;a1=a0+1;N=0;while abs(a1-a0)>ea=a0-subs(f,a0)/subs(df,a0); a1=a0; a0=a; N=N+1; endfprintf('a=%0.6f',a) N运行结果: >> newtonplease enter your equation:exp(x)+2^(-x)+2*cos(x)-6 f =exp(x)+2^(-x)+2*cos(x)-6 please enter you x(0):2df =exp(x)-2^(-x)*log(2)-2*sin(x) a=1.829384 N =42.利用Newton 迭代法求多项式07951.2954.856.104.5x 234=+-+-x x x的所有实零点,注意重根的问题。
数值分析2024上机实验报告数值分析是计算数学的一个重要分支,它研究如何用数值方法来解决数学问题。
在数值分析的学习过程中,学生需要通过上机实验来巩固理论知识,并学会使用相应的数值方法来解决实际问题。
本篇报告将详细介绍2024年度数值分析上机实验的内容和结果。
一、实验内容2024年度数值分析上机实验分为四个部分,分别是:方程求根、插值与拟合、数值积分和常微分方程的数值解。
1.方程求根这部分实验要求使用数值方法求解给定的非线性方程的根。
常见的数值方法有二分法、牛顿法、割线法等。
在实验过程中,我们需要熟悉这些数值方法的原理和实现步骤,并对不同方法的收敛性进行分析和比较。
2.插值与拟合这部分实验要求使用插值和拟合方法对给定的一组数据进行拟合。
插值方法包括拉格朗日插值、牛顿插值等;拟合方法包括最小二乘拟合、多项式拟合等。
在实验中,我们需要熟悉插值和拟合方法的原理和实现步骤,并对不同方法的精度和稳定性进行比较。
3.数值积分这部分实验要求使用数值方法计算给定函数的积分。
常见的数值积分方法有梯形法则、辛普森法则、龙贝格积分等。
在实验过程中,我们需要熟悉这些数值积分方法的原理和实现步骤,并对不同方法的精度和效率进行比较。
4.常微分方程的数值解这部分实验要求使用数值方法求解给定的常微分方程初值问题。
常见的数值方法有欧拉法、改进的欧拉法、四阶龙格-库塔法等。
在实验中,我们需要熟悉这些数值解方法的原理和实现步骤,并对不同方法的精度和稳定性进行比较。
二、实验结果在完成2024年度数值分析上机实验后,我们得到了以下实验结果:1.方程求根我们实现了二分法、牛顿法和割线法,并对比了它们的收敛速度和稳定性。
结果表明,割线法的收敛速度最快,但在一些情况下可能会出现振荡;二分法和牛顿法的收敛速度相对较慢,但稳定性较好。
2.插值与拟合我们实现了拉格朗日插值和最小二乘拟合,并对比了它们的拟合效果和精度。
结果表明,拉格朗日插值在小区间上拟合效果较好,但在大区间上可能出现振荡;最小二乘拟合在整体上拟合效果较好,但可能出现过拟合。
第一题:1、已知A 与b12.38412 2.115237 -1.061074 1.112336 -0.1135840.718719 1.742382 3.067813 -2.031743 2.11523719.141823 -3.125432 -1.012345 2.189736 1.563849-0.784165 1.112348 3.123124 -1.061074 -3.125A =43215.567914 3.123848 2.031454 1.836742-1.056781 0.336993 -1.010103 1.112336 -1.012345 3.12384827.108437 4.101011-3.741856 2.101023 -0.71828 -0.037585 -0.1135842.189736 2.031454 4.10101119.8979180.431637-3.111223 2.121314 1.784317 0.718719 1.563849 1.836742 -3.741856 0.4316379.789365-0.103458 -1.103456 0.238417 1.742382 -0.784165 -1.056781 2.101023-3.111223-0.10345814.7138465 3.123789 -2.213474 3.067813 1.112348 0.336993-0.71828 2.121314-1.103456 3.12378930.719334 4.446782 -2.031743 3.123124 -1.010103-0.037585 1.7843170.238417-2.213474 4.44678240.00001[ 2.1874369 33.992318 -25.173417 0.84671695 1.784317 -86.612343 1.1101230 4.719345 -5.6784392]TB ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦=(1)用Househloser 变换,把A 化为三对角阵(并打印B )。
数值分析第一次上机练习实验报告一、实验目的本次实验旨在通过上机练习,加深对数值分析方法的理解,并掌握实际应用中的数值计算方法。
二、实验内容1. 数值计算的基本概念和方法在本次实验中,我们首先回顾了数值计算的基本概念和方法。
数值计算是一种通过计算机进行数值近似的方法,其包括近似解的计算、误差分析和稳定性分析等内容。
2. 方程求解的数值方法接下来,我们学习了方程求解的数值方法。
方程求解是数值分析中非常重要的一部分,其目的是找到方程的实数或复数解。
我们学习了二分法、牛顿法和割线法等常用的数值求解方法,并对它们的原理和步骤进行了理论学习。
3. 插值和拟合插值和拟合是数值分析中常用的数值逼近方法。
在本次实验中,我们学习了插值和拟合的基本原理,并介绍了常见的插值方法,如拉格朗日插值和牛顿插值。
我们还学习了最小二乘拟合方法,如线性拟合和多项式拟合方法。
4. 数值积分和数值微分数值积分和数值微分是数值分析中的两个重要内容。
在本次实验中,我们学习了数值积分和数值微分的基本原理,并介绍了常用的数值积分方法,如梯形法和辛卜生公式。
我们还学习了数值微分的数值方法,如差商法和牛顿插值法。
5. 常微分方程的数值解法常微分方程是物理和工程问题中常见的数学模型,在本次实验中,我们学习了常微分方程的数值解法,包括欧拉法和四阶龙格-库塔法。
我们学习了这些方法的步骤和原理,并通过具体的实例进行了演示。
三、实验结果及分析通过本次实验,我们深入理解了数值分析的基本原理和方法。
我们通过实际操作,掌握了方程求解、插值和拟合、数值积分和数值微分以及常微分方程的数值解法等数值计算方法。
实验结果表明,在使用数值计算方法时,我们要注意误差的控制和结果的稳定性。
根据实验结果,我们可以对计算结果进行误差分析,并选择适当的数值方法和参数来提高计算的精度和稳定性。
此外,在实际应用中,我们还需要根据具体问题的特点和条件选择合适的数值方法和算法。
四、实验总结通过本次实验,我们对数值分析的基本原理和方法有了更加深入的了解。
#### 一、实验目的本次实验旨在通过MATLAB软件,对矩阵进行数值计算,掌握矩阵的基本操作、运算函数的使用,以及解决实际问题的能力。
通过实验,加深对线性代数基本理论的理解,提高数值计算技能。
#### 二、实验环境软件:MATLAB R2020a硬件:****************************,8GB RAM#### 三、实验内容1. 矩阵的创建与操作(1)创建矩阵:通过MATLAB内置函数创建不同类型的矩阵,如`zeros`、`ones`、`rand`等。
```matlabA = zeros(3,3); % 创建3x3零矩阵B = ones(2,2); % 创建2x2单位矩阵C = rand(4,4); % 创建4x4随机矩阵```(2)矩阵的引用:使用矩阵的行和列索引访问矩阵元素。
```matlabE = A(1,1); % 访问矩阵A的第一个元素```(3)矩阵的运算:进行矩阵的加法、减法、乘法等运算。
```matlabD = A + B; % 矩阵A和B相加F = A . C; % 矩阵A和C对应元素相乘```2. 矩阵的基本运算(1)矩阵的逆:计算矩阵的逆矩阵。
```matlabA_inv = inv(A);```(2)矩阵的行列式:计算矩阵的行列式值。
```matlabdet_A = det(A);```(3)矩阵的秩:计算矩阵的秩。
```matlabrank_A = rank(A);```(4)矩阵的迹:计算矩阵的迹。
```matlabtrace_A = trace(A);```3. 矩阵分解(1)奇异值分解(SVD):对矩阵进行奇异值分解。
```matlab[U, S, V] = svd(A);```(2)LU分解:将矩阵分解为下三角矩阵和上三角矩阵。
```matlab[L, U] = lu(A);```4. 解线性方程组(1)使用矩阵的逆解方程组。
```matlabb = [1; 2; 3];x = A_inv b;```(2)使用矩阵分解方法解方程组。
《矩阵与数值分析》上机大作业1.给定n 阶方程组Ax b =,其中6186186186A ⎛⎫ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭ ,7151514b ⎛⎫ ⎪ ⎪ ⎪= ⎪⎪⎪⎝⎭则方程组有解(1,1,,1)Tx = 。
对10n =和84n =,分别用Gauss 消去法和列主元消去法解方程组,并比较计算结果。
%产生三对角矩阵 n=84; %或n=10;A=zeros(n); b=zeros(1,n); for i=1:n-1A(i,i)=6;A(i,i+1)=1;A(i+1,i)=8; endA(n,n)=6;for i=2:n-1 b(1)=6; b(i)=15; b(n)=14; end Ab=[A b'];%Gauss 消元法for j=1:n-1 %按列循环 for k=j+1:n %消元Ab(k,:)=Ab(k,:)-Ab(j,:)*(Ab(k,j)/Ab(j,j)); end endx(n)=Ab(n,n+1)/Ab(n,n); for i=n-1:-1:1 %回代法求x for j=n:-1:i+1Ab(i,n+1)=Ab(i,n+1)-Ab(i,j)*x(j); endx(i)=Ab(i,n+1)/Ab(i,i); end(1)当n=10时,Gauss 消去法 Gauss 列主元消去法 x=1.000000000000000 x=1.000000000000000 1.000000000000000 1.0000000000000001.000000000000000 1.000000000000000 1.000000000000001 1.000000000000000 0.999999999999998 1.000000000000000 1.000000000000004 1.000000000000000 0.999999999999993 1.000000000000000 1.000000000000012 1.000000000000000 0.999999999999979 1.000000000000000 1.000000000000028 1.000000000000000(2) 当n=84时,Gauss 消去法的解是错解Columns 34 through 392147483649.00000 -4294967295.00000 8589934592.99999 -17179869182.9999 34359738368.9998Gauss 列主元消去法x 与x=(1,1…1)T 偏差不大 Columns 34 through 391.000000172108412 0.999999661246936 1.000000655651093 0.999998776117961 1.000002098083496综上,高斯列主元消去法可以避免小数作除数带来的误差,获得满意的数值解。
矩阵数值分析上机实习报告专业:信息与计算科学班级:学号:姓名:目录一、Doolittle分解 (1)二、Jacobi算法 (5)三、逆幂法 (8)四、古典雅克比 (12)五、割线法 (16)一、Doolittle 分解1、目的意义:在Gauss 消元法中,对于n 阶方程组,应用消去发经过n-1步消元之后,得到一个与Ax=b 等价的代数线性方程组)1()1(--=n n b x A ,而且)1(-n A 为一个上三角矩阵.所以我们想是否能把矩阵A 分解成一个下三角阵与一个上三角阵的乘积 A=LR,其中L 为下三角阵,R 为上三角阵.就变成了两个三角形方程组⎩⎨⎧==y Rx b Ly ,的求解问题.2、算法:Setp1:利用for 循环求出11[][][][][][]k p r k j a k j l kp r kp -==-∑,11[]([][][])/[][]k p l ik a ik l ip r ip r k k -==-∑Step2:11[][][][]i k y i b i l ik y k -==-∑得出1[]([][][])/[][]nk i x i y i r ik x k r i i =+=-∑3、流程图4、源程序:#include"stdio.h"#define N 10void main(){float a[N][N],b[N],x[N],y[N],L[N][N],r[N][N],sum1,sum2,sum3,sum4;int i,j,k,p,n;printf("*****Doolittle分解法求方程组*****\n");printf("输入n:");scanf("%d",&n);printf("请输入系数矩阵:\n");for(i=1;i<=n;i++){for(j=1;j<=n;j++)scanf("%f",&a[i][j]);}printf("请输入等号右边的值:\n");for(i=1;i<=n;i++)scanf("%f",&b[i]);for(i=1;i<=n;i++) //给系数矩阵赋初值{for(j=1;j<=n;j++){if(i==j){L[i][j]=1;r[i][j]=0;}else{L[i][j]=0;r[i][j]=0;}}}for (k=1;k<=n;k++){for(i=k;i<=n;i++){sum1=0;for(p=1;p<=k-1;p++){sum1+=L[k][p]*r[p][i];}r[k][i]=a[k][i]-sum1;}for(i=k+1;i<=n;i++){sum2=0;for(p=1;p<=k-1;p++){sum2+=L[i][p]*r[p][k];}L[i][k]=(a[i][k]-sum2)/r[k][k];}}printf("输出L矩阵:\n");for(i=1;i<=n;i++){for(j=1;j<=n;j++)printf("%f\t",L[i][j]);printf("\n");}printf("输出R矩阵:\n");for(i=1;i<=n;i++){for(j=1;j<=n;j++)printf("%f\t",r[i][j]);printf("\n");}printf("输出y\n");for(i=1;i<=n;i++){sum3=0;for(j=1;j<i;j++)sum3=sum3+L[i][j]*y[j];y[i]=b[i]-sum3;}for(i=1;i<=n;i++)printf("y[%d]=%f\n",i,y[i]);printf("\n");printf("输出x\n");for(i=n;i>=1;i--){sum4=0;for(k=i+1;k<=n;k++)sum4=sum4+r[i][k]*x[k];x[i]=(y[i]-sum4)/r[i][i];printf("x[%d]=%f\t",i,x[i]);}}5、程序输出:*****Doolittle分解法求方程组*****输入n:3请输入系数矩阵:2 1 11 3 21 2 2请输入等号右边的值:4 6 5输出L矩阵:1.000000 0.000000 0.000000 0.500000 1.000000 0.000000 0.500000 0.600000 1.000000 输出R矩阵:2.000000 1.000000 1.0000000.000000 2.500000 1.5000000.000000 0.000000 0.600000输出yy[1]=4.000000y[2]=4.000000y[3]=0.600000输出xx[3]=1.000000 x[2]=1.000000 x[1]=1.000000 Press any key to continue6、参考文献:[1]刑志栋,矩阵数值分析,陕西:陕西科学技术出版社, 2005。
二、Jacobi 算法1、问题描述设方程组Ax=b 的系数矩阵A 非奇异而且),...,2,1(0n i a ii =≠,将A 分裂为A=D+L+U,可以使计算简便。
其中),,...,,(2211nn a a a diag D =⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=0...............0...00 (00)2121n n a a a L ,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=0...00 00...02112n n a a a U2、思想:A=D+L+U ,其中),,...,,(2211nn a a a diag D =⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=0...............0...00 (002)121n n a a a L ,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=0...00 00...02112n n a a a U将方程组n ,...,2,1i ,b x a i n1j j ij ==∑=乘以ii a 1,得到等价的方程组⎪⎪⎪⎭⎫⎝⎛-=∑≠=n i j 1j j ij i ii i x a b a 1x ,i=1,2,…n3、算法:Step1:给定一组x ,即初值。
Step2:用for 循环计算: x[k+1]=(b[i]-∑∑+=-=-n1i j 1i 1j ]j [x ]j ][i [a ]j [x ]j ][i [a )/a[k][k].Step3:当abs(x[k+1]-x[k])<eps 时停止。
4、程 序: #include"stdio.h" #include"math.h" #define N 3#define M 20#define EPS 1e-6void main(){floata[N][N],b[N],x[M][N],l,p;int i,j,k;printf("请输入系数矩阵:\n");for(i=1;i<=N;i++){for(j=1;j<=N;j++){scanf("%f",&a[i][j]);}}printf("请输入方程组等号右边的值:\n");for(i=1;i<=N;i++)scanf("%f",&b[i]);printf("请输入初始变量:\n");for(i=1;i<=N;i++)scanf("%f",&x[0][i]);for(k=0;k<=M;k++){p=0;printf("第%d次迭代:\n",k+1);for(i=1;i<=N;i++){l=0;for(j=1;j<=N;j++){if(j==i) l=l;elsel+=a[i][j]*x[k][j];}x[k+1][i]=(b[i]-l)/a[i][i];p+=fabs(x[k+1][i]-x[k][i]);printf("x[%d]=%f\n",i,x[k+1][i]);}if(p<=EPS)break;}}5、算法输出:请输入系数矩阵:1 2 -21 1 12 2 1请输入方程组等号右边的值: 1 3 5请输入初始变量:0 0 0第1次迭代:x[1]=1.000000x[2]=3.000000x[3]=5.000000第2次迭代:x[1]=5.000000x[2]=-3.000000x[3]=-3.000000第3次迭代:x[1]=1.000000x[2]=1.000000x[3]=1.000000第4次迭代:x[1]=1.000000x[2]=1.000000x[3]=1.000000Press any key to continue6、参考文献[1]刑志栋,矩阵数值分析,陕西:陕西科学技术出版社, 2005。
三、逆幂法1、问题描述逆幂法用于计算奇异矩阵A 的按模最小特征值和相应的特征向量。
2、思想step1:1-=k k z y L step2: k k y Uy =step3: i k ni k y m )(max 1≤≤= k=1,2,…step4: k k k m y z /=其中LU=A.L 是单位下三角矩阵,U 为上三角矩阵。
3、流程图4、源程序: #include<stdio.h> #define exp 0.001 #define N 10void main(){int i,j,k,n,v;double a[N+1][N+1],c[N+1][N+1],d[N+1][N+1],l[N+1][N+1],r[N+1][N+1];double x[N+1],y[N+1],z[N+1],b[N+1],h[N+1],g[N+1],m[N+1];double p,q,t,s,u1,u2,w,max;printf("请输入n:");scanf("%d",&n);printf("请输入特征值u1:");scanf("%lf",&u1);printf("请输入初始向量z[]:\n");for(i=1;i<=n;i++){scanf("%lf",&z[i]); }printf("请输入系数矩阵c[][]:\n"); for(i=1;i<=n;i++)for(j=1;j<=n;j++){scanf("%lf",&c[i][j]);if(i==j)d[i][j]=1;else d[i][j]=0;a[i][j]=c[i][j]-u1*d[i][j];}for(i=1;i<=n;i++){for(j=i;j<=n;j++){for(k=1,p=0;k<=i-1;k++)p=p+l[i][k]*r[k][j];r[i][j]=a[i][j]-p;}for(j=i+1;j<=n;j++){for(k=1,q=0;k<=i-1;k++)q=q+l[j][k]*r[k][i];l[j][i]=(a[i][j]-q)/r[i][i];}l[i][i]=1;}for(i=1;i<=n;i++){m[i]=0;for(j=i;j<=n;j++)printf("r[%d][%d]:%lf\n",i,j,r[i][j]);for(j=i+1;j<=n;j++) printf("l[%d][%d]:%lf\n",j,i,l[j][i]); }for(v=0; ;v++){printf("\n",v);w=0;max=0;for(i=1;i<=n;i++){m[i]=x[i];for(k=1,t=0;k<=i-1;k++)t=t+l[i][k]*y[k];y[i]=z[i]-t;}printf("y[]=");for(i=1;i<=n;i++)printf("%lf\t",y[i]);printf("\n");printf("z[]=");for(i=1;i<=n;i++)printf("%lf\t",z[i]);printf("\n");printf("x[]=");for(i=n;i>=1;i--){for(k=i+1,s=0;k<=n;k++)s=s+r[i][k]*x[k];x[i]=(y[i]-s)/r[i][i];printf("%lf\t",x[i]);if(x[i]>=0)g[i]=x[i];if(x[i]<0)g[i]=-x[i];if(max<g[i])max=g[i];}printf("\n");for(i=1;i<=n;i++) {z[i]=x[i]/max;b[i]=m[i]-x[i];if(b[i]>=0)h[i]=b[i];else h[i]=-b[i];w=w+h[i];}printf("z[]=");for(i=1;i<=n;i++)printf("%lf\t",z[i]);printf("\n");if(w<exp)break;}u2=1/max+u1;printf("u2=%lf\n",u2);}5、运行结果:6、参考文献[1]刑志栋,矩阵数值分析,陕西:陕西科学技术出版社, 2005。