数值分析实验报告(Matlab实现)
学 生 实 验 报 告 实验课程名称 数值分析 开课实验室 数学与统计学院实验室 学 院 2010 年级 数学与应用数学专业班 01班 学 生 姓 名 学 号 开 课 时 间 2012 至 2013 学年第 一 学期
总 成 绩 教师签名 1
课程 名称 数值分析 实验项目 名 称 Gauss消元法
实验项目类型
验证 演示 综合 设计 其他 指导 教师 何光辉 成 绩 √
一、 实验目的: (1)高斯列主元消去法求解线性方程组的过程 (2)熟悉用迭代法求解线性方程组的过程 (3)设计出相应的算法,编制相应的函数子程序
二、 实验内容 分别用高斯列主元消元法和直接消元法求解线性方程组:
72510139144432113124330102
4321
xxx
x
三、 实验原理 对于线性方程组
1111221121122222n1122n=b=b=bnnnn
nnnn
axaxaxaxaxaxaxaxax……………
(1)
常记为矩阵形式 Axb (2)
根据高等代数的知识,若0A,上式的解存在且唯一。 (1) Gauss直接消元法 考虑上述线性方程组的增广矩阵[:]Ab,对增广矩阵进行行变换,将(2)式化为等价的三角形方阵,然后回代解之,这就是Gauss消元法。具体如下: a) 消元
①令(1),,1,2nijijaaij,…,;(1),1,2niibbi,…,
②对k=1到n-1,若()0kkka,进行 2
()()(1)(1)()()(1)(),1,2,n0,1,2,n,,1,2,n,1,2,nkikikk
kk
kikkkkijijikkjkkkiiikkalikkaaikkaalaijkkbblbikk
…,
…,…,…,
b) 回代,若a0nnn
()()11()nnnnnnniiiiijji
jkiibxbxbaxa
(2) Gauss列主元消元法 设列主元消元法已完成Axb的第k-1(11kn)次消元,的到方程组
()()kkAxbAxb
在进行第k次消元前,先进行2个步骤: a) 在()kkka至()knka这一列内选出最大值,即()(),maxkkkikikkinaa,若(),0kkika,此时0A
方程组无确定解,应给出退出信息。 b) 若(),0kkika,则交换第ki行和k行,然后用Gauss消元法进行消元。
四、 MATLAB软件实现 (1) 写出Gauss消元法和列主元消元法实现的MATLAB函数 根据以上的算法,写出如下程序: %%%%%%%Gauss消元法%%%%%%%%%%%% function y=Gauss1(A,b) [m,n]=size(A); %检查系数正确性 if m~=n error('矩阵A的行数和列数必须相同'); return; end if m~=size(b) error(' b的大小必须和A的行数或A的列数相同'); return; end %再检查方程是否存在唯一解 if rank(A)~=rank([A,b]) error(' A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解'); return; 3
end %这里采用增广矩阵行变换的方式求解 c=n+1; A(:,c)=b; %%消元过程 for k=1:n-1 A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); end %%回代结果 x=zeros(length(b),1); x(n)=A(n,c)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k); end %显示计算结果 %disp('x='); %disp(x); y=x;
% %%%%%%%%%%%高斯列主元消元法求解线性方程组Ax=b%%%%%%%%%%%%%% %A为输入矩阵系数,b为方程组右端系数 %方程组的解保存在x变量中 function y=Gauss_line(A,b) format long;% 设置为长格式显示,显示15位小数 [m,n]=size(A); %先检查系数正确性 if m~=n error('矩阵A的行数和列数必须相同'); return; end if m~=size(b) error(' b的大小必须和A的行数或A的列数相同'); return; end %再检查方程是否存在唯一解 if rank(A)~=rank([A,b]) error(' A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解'); return; end c=n+1; A(:,c)=b; %(增广) for k=1:n-1 [r,m]=max(abs(A(k:n,k))); %选主元 m=m+k-1; %修正操作行的值 4
if(A(m,k)~=0) if(m~=k) A([k m],:)=A([m k],:); %换行 end A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); %消去 end end x=zeros(length(b),1); %回代求解 x(n)=A(n,c)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k); end y=x; format short;% 设置为默认格式显示,显示5位
(2) 建立MATLAB界面 利用MATLAB的GUI建立如下界面求解线性方程组:
详见程序。 五、 计算实例、数据、结果、分析 下面我们对以上的结果进行测试,求解: 5
72510139144432113124330102
4321
xxx
x
输入数据后点击和,得到如下结果:
更改以上数据进行测试,求解如下方程组: 1234
43211343212343112341xxxx
得到如下结果: 6
六、 实验中遇到的问题及解决办法 在本实验中,遇到的问题主要有两个: (1) 如何将上述的Gauss消元法的算法在MATLAB中实现 针对此问题我借鉴了网上以及 课本上的算法的MATLAB实现的程序; (2) 如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解 针对此问题,我通过网上的一些关于MATLAB的GUI设计的相关资料,总结经验完成了此项任务。
七、 实验结论 通过以上的测试,我们发现以上算法和程序能够求出线性方程组的比较精确解。
八、 参考文献 [1]杨大地,王开荣 .2006.数值分析.北京:科学出版社 [2]何光辉.2008. 数值分析实验. 重庆大学数理学院数学实验教学中心 [3]百度文库,百度知道
教师签名 年 月 日
课程 名称 数值分析 实验项目 名 称 插值方法
实验项目类型
验证 演示 综合 设计 其他 指导 教师 何光辉 成 绩 √
一、 实验目的: (1) 学会拉格朗日插值、牛顿插值等基本方法 (2) 设计出相应的算法,编制相应的函数子程序 (3) 会用这些函数解决实际问题
二、 实验内容 (1)设计拉格朗日插值算法,编制并调试相应的函数子程序 (2)设计牛顿插值算法,编制并调试相应的函数子程序 (3)给定函数四个点的数据如下: X 1.1 2.3 3.9 5.1 Y 3.887 4.276 4.651 2.117 试用拉格朗日插值确定函数在x=2.101,4.234处的函数值。
(4)已知,,,392411用牛顿插值公式求5的近似值。
三、 实验原理 (1) 拉格朗日插值
数值分析实验报告2
实验报告实验项目名称函数逼近与快速傅里叶变换实验室数学实验室所属课程名称数值逼近实验类型算法设计实验日期班级学号姓名成绩512*x^10 - 1280*x^8 + 1120*x^6 - 400*x^4 + 50*x^2 - 1并得到Figure,图像如下:实验二:编写程序实现[-1,1]上n阶勒让德多项式,并作画(n=0,1,…,10 在一个figure中)。
要求:输入Legendre(-1,1,n),输出如a n x n+a n-1x n-1+…多项式。
在MATLAB的Editor中建立一个M-文件,输入程序代码,实现勒让德多项式的程序代码如下:function Pn=Legendre(n,x)syms x;if n==0Pn=1;else if n==1Pn=x;else Pn=expand((2*n-1)*x*Legendre(n-1)-(n-1)*Legendre(n-2))/(n);endx=[-1:0.1:1];A=sym2poly(Pn);yn=polyval(A,x);plot (x,yn,'-o');hold onend在command Windows中输入命令:Legendre(10),得出的结果为:Legendre(10)ans =(46189*x^10)/256 - (109395*x^8)/256 + (45045*x^6)/128 - (15015*x^4)/128 + (3465*x^2)/256 - 63/256并得到Figure,图像如下:实验三:利用切比雪夫零点做拉格朗日插值,并与以前拉格朗日插值结果比较。
在MATLAB的Editor中建立一个M-文件,输入程序代码,实现拉格朗日插值多项式的程序代码如下:function [C,D]=lagr1(X,Y)n=length(X);D=zeros(n,n);D(:,1)=Y';for j=2:nfor k=j:nD(k,j)=(D(k,j-1)- D(k-1,j-1))/(X(k)-X(k-j+1));endendC=D(n,n);for k=(n-1):-1:1C=conv(C,poly(X(k)));m=length(C);C(m)= C(m)+D(k,k);end在command Windows 中输入如下命令:clear,clf,hold on;k=0:10;X=cos(((21-2*k)*pi)./22); %这是切比雪夫的零点Y=1./(1+25*X.^2);[C,D]=lagr1(X,Y);x=-1:0.01:1;y=polyval(C,x);plot(x,y,X,Y,'.');grid on;xp=-1:0.01:1;z=1./(1+25*xp.^2);plot(xp,z,'r')得到Figure ,图像如下所示:比较后发现,使用切比雪夫零点做拉格朗日插值不会发生龙格现象。
数值分析实验报告
实验2.1 多项式插值的振荡现象实验目的:在一个固定的区间上用插值逼近一个函数,显然Lagrange 插值中使用的节点越多,插值多项式的次数就越高。
我们自然关心插值多项式的次数增加时,Ln(x)是否也更加靠近被逼近的函数。
Runge 给出的一个例子是极著名并富有启发性的。
实验容:设区间[-1,1]上函数 f(x)=1/(1+25x 2)。
考虑区间[-1,1]的一个等距划分,分点为 x i = -1 + 2i/n ,i=0,1,2,…,n ,则拉格朗日插值多项式为201()()125nn i i i L x l x x ==+∑. 其中,l i (x),i=0,1,2,…,n 是n 次Lagrange 插值基函数。
实验步骤与结果分析:实验源程序function Chap2Interpolation% 数值实验二:“实验2.1:多项式插值的震荡现象”% 输入:函数式选择,插值结点数% 输出:拟合函数及原函数的图形promps = {'请选择实验函数,若选f(x),请输入f,若选h(x),请输入h,若选g(x),请输入g:'};titles = 'charpt_2';result = inputdlg(promps,'charpt 2',1,{'f'});Nb_f = char(result);if(Nb_f ~= 'f' & Nb_f ~= 'h' & Nb_f ~= 'g')errordlg('实验函数选择错误!');return;endresult = inputdlg({'请输入插值结点数N:'},'charpt_2',1,{'10'});Nd = str2num(char(result));if(Nd <1)errordlg('结点输入错误!');return;endswitch Nb_fcase 'f'f=inline('1./(1+25*x.^2)'); a = -1;b = 1;case 'h'f=inline('x./(1+x.^4)'); a = -5; b = 5;case 'g'f=inline('atan(x)'); a = -5; b= 5;endx0 = linspace(a, b, Nd+1); y0 = feval(f, x0);x = a:0.1:b; y = Lagrange(x0, y0, x);fplot(f, [a b], 'co');hold on;plot(x, y, 'b--');xlabel('x'); ylabel('y = f(x) o and y = Ln(x)--');%--------------------------------------------------------------------function y=Lagrange(x0, y0, x);n= length(x0); m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif(j ~= k)p = p*(z - x0(j))/(x0(k) - x0(j));endends = s + p*y0(k);endy(i) = s;end实验结果分析(1)增大分点n=2,3,…时,拉格朗日插值函数曲线如图所示。
数值分析实验报告--解线性方程组的迭代法及其并行算法
disp('请注意:高斯-塞德尔迭代的结果没有达 到给定的精度,并且迭代次数已经超过最大迭 代次数max1,方程组的精确解jX和迭代向量X 如下: ') X=X';jX=jX' end end X=X';D,U,L,jX=jX'
高斯-塞德尔的输入为:
A=[10 2 3;2 10 1;3 1 10]; b=[1;1;2]; X0=[0 0 0]'; X=gsdddy(A,b,X0,inf, 0.001,100) A=[10 2 3;2 10 1;3 1 10]; 请注意:因为对角矩阵 D 非奇异,所以此方程组有解.
0.0301 0.0758 0.1834
8.心得体会:
这已经是第三次实验了, 或多或少我已经对 MATLAB 有了更多的了 解与深入的学习。通过这次实验我了解了雅可比迭代法和高斯- 塞德尔迭代法的基本思想,虽然我们不能熟练编出程序,但还是 能看明白的。运行起来也比较容易,让我跟好的了解迭代法的多 样性,使平常手算的题能得到很好的验证。通过这次实验让我对 MATLAB 又有了更深一层的认识,使我对这门课兴趣也更加浓厚。
运行雅可比迭代程序输入: A=[10
b=[1;1;2];X0=[0 0 0]'; X=jacdd(A,b,X0,inf,0.001,100)
2 3;2 10 1;3 1 10];
结果为:
k= 1 X=
0.1000 k= 2 X= 0.0200 k= 3 X= 0.0400 k= 4 X= 0.0276 k= 5 X= 0.0314 k= 6 X= 0.0294 k= 7 X= 0.0301 k= 8 X= 0.0297
6、 设计思想:先化简,把对角线的项提到左边,其它项
数值分析实验2014
数值分析实验(2014,9,16~10,28)信计1201班,人数34人数学系机房数值分析计算实习报告册专业__________________学号_______________姓名_______________2014~2015年第一学期实验一数值计算的工具Matlab1. 解释下MATLABS序的输出结果程序:t=0.1n=1:10e=n/10-n*te 的结果:0 0 -5.5511e-017 0 0-1.1102e-016 -1.1102e-016 0 0 02. 下面MATLABS序的的功能是什么?程序:x=1;while 1+x>1,x=x/2,pause(0.02),e nd用迭代法求出x=x/2,的最小值x=1;while x+x>x,x=2*x,pause(0.02),e nd用迭代法求出x=2*x,的值,使得2x>Xx=1;while x+x>x,x=x/2,pause(0.02),e nd用迭代法求出x=x/2,的最小值,使得2x>X3. 考虑下面二次代数方程的求解问题2ax bx c = 0公式x=电上4ac是熟知的,与之等价地有_____________________________ ,对于2a-b ■ b -4aca =1,b =100000000,c =1,应当如何选择算法。
b ~4ac计算,因为b与b2— 4ac相近,两个相加减不宜应该用2a u做分母3 5 74. 函数sin(x)有幂级数展开sin x = x - x - - ■■3! 5! 7!利用幕级数计算sinx的MATLAB程序为fun cti on s=powers in(x)s=0;t=x;n=1;while s+t~=s;s=s+t ;t=-x A2/ ((n+1)*(n+2) ) *t ;n=n+2 ;endt仁cputime;pause(10);t2=cputime;t0=t2-t1(a) 解释上述程序的终止准则。
数值分析Hilbert矩阵病态线性方程组的求解Matlab程序
(Hilbert矩阵)病态线性方程组的求解理论分析表明,数值求解病态线性方程组很困难。
考虑求解如下的线性方程组的求解1Hx=b,期中H是Hilbert矩阵,H=(h j)四,h j=,।,j=1,2,…,nij-11,估计矩阵的2条件数和阶数的关系2 .对不同的n,取x=(1,1,…,1)w n,分别用Gauss消去,Jacobi迭代,Gauss-seidel迭代,SOR迭代和共轲梯度法求解,比较结果。
3 .结合计算结果,试讨论病态线性方程组的求解。
第1小题:condition.m。
蟆1小题程序t1=20;%>数n=20x1=1:t1;y1=1:t1;fori=1:t1H=hilb(i);y1(i)=log(cond(H));endplot(x1,y1);xlabel('阶数n');ylabel('2-条件数的对数(log(cond(H))');title('2-条件数的对数(10g(cond(H))与阶数n的关系图’);t2=200;%J>数n=200x2=1:t2;y2=1:t2;fori=1:t2H=hilb(i);y2(i)=log(cond(H));endplot(x2,y2);xlabel('阶数n');ylabel('2-条件数的对数(log(cond(H))');title('2-条件数的对数(10g(cond(H))与阶数n的关系图’);画出Hilbert矩阵2-条件数的对数和阶数的关系n=200时n=20时马KiMn的美方四*--wficEJ-帛=«』从图中可以看出,1)在n小于等于13之前,图像近似直线log(cond(H))-1.519n-1.8332)在n大于13之后,图像趋于平缓,并在一定范围内上下波动,同时随着n的增加稍有上升的趋势第2小题:soke.m%蹴2小题主程序N=4000;xGauss=zeros(N,1);xJacobi=zeros(N,1);xnJ=zeros(N,1);xGS=zeros(N,1);xnGS=zeros(N,1);xSOR=zeros(N,1);xnSOR=zeros(N,1);xCG=zeros(N,1);xnCG=zeros(N,1);forn=1:N;x=ones(n,1);t=1.1;%初始值偏差x0=t*x;%迭代初始值e=1.0e-8;%^定的误差A=hilb(n);b=A*x;max=100000000000;%可能最大的迭代次数w=0.5;%SORt代的松弛因子G=Gauss(A,b);[J,nJ]=Jacobi(A,b,x0,e,max);[GS,nGS]=G_S(A,b,x0,e,max);[S_R,nS_R]=SOR(A,b,x0,e,max,w);[C_G,nC_G]=CG(A,b,x0,e,max);normG=norm(G'-x);xGauss(n)=normG;normJ=norm(J-x);nJ;xJacobi(n)=normJ;xnJ(n)=nJ;normGS=norm(GS-x);nGS;xGS(n)=normGS;xnGS(n)=nGS;normS_R=norm(S_R-x);nS_R;xSOR(n)=normS_R;xnSOR(n)=nS_R;normC_G=norm(C_G-x);nC_G;xCG(n)=normC_G;xnCG(n)=nC_G;endGauss.m%Gauss消去法functionx=Gauss(A,b)n=length(b);l=zeros(n,n);x=zeros(1,n);脸肖去过程fori=1:n-1forj=i+1:nl(j,i)=A(j,i)/A(i,i);fork=i:nA(j,k)=A(j,k)-l(j,i)*A(i,k);endb(j)=b(j)-l(j,i)*b(i);endend%回代过程x(n)=b(n)/A(n,n);fori=n-1:-1:1c=A(i,:).*x;x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);endJacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function[x,n]=Jacobi(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;whilenorm(x-x0)>ex0=x;x=B*x0+f;n=n+1;ifn>mdisp('Jacobi迭代次数过多,迭代可能不收敛‘);break;endendG_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function[x,n]=G_S(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;whilenorm(x-x0)>ex0=x;x=B*x0+f;n=n+1;ifn>mdisp('Gauss-Seidel迭代次数过多,迭代可能不收敛’);break;endendSOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子function[x,n]=SOR(A,b,x0,e,m,w)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=B*x0+f;n=1;whilenorm(x-x0)>ex0=x;x=B*x0+f;n=n+1;ifn>mdisp('SOR超松弛迭代次数过多,迭代可能不收敛’);break;endendCG.m%C献朝梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function[x,n]=CG(A,b,x0,e,m)r=b-A*x0;p=r;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=1;whilenorm(r1)>ebelta=(r1'*r1)/(r'*r);p=r1+belta*p;r=r1;x0=x;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=n+1;ifn>mdisp('CG共朝梯度法迭代次数过多,迭代可能不收敛’);break;endend第2小题选择不同的阶数n,取x=(l,l,…RI分别使用Gauss消去,Jacobi迭代,Gauss-seidel迭代,SOR迭代和共挽梯度法(CG)求解,比较结果。
MATLAB软件及高斯勒让德求积公式
0.5689 0.2369 0.4786 0.2369 0.4786
x =
0 0.9062 0.5385 -0.9062 -0.5385
>> [A,x]=Guass1(5)
% a,b·Ö±ðÊÇ»ý·ÖµÄÉÏÏÂÏÞ£»
% n+1Ϊ½Úµã¸öÊý£»
% mÊǵ÷ÓÃf1.mÖеڼ¸¸ö±»»ýº¯Êý£»
[A,x]=Guass1(n)
g=0;
fori=1:n+1
y(i)=(b-a)/2*x(i)+(a+b)/2;
f(i)=f1(m,y(i));
g=g+(b-a)/2*f(i)*A(i);
f=diff(f,i);
t=solve(f);
forj=1:i
fork=1:i
X(j,k)=t(k)^(j-1);
end
ifmod(j,2)==0
B(j)=0;
else
B(j)=2/j;
end
end
X=inv(X);
forj=1:i
A(j)=0;
x(j)=0;
fork=1:i
A(j)=A(j)+X(j,k)*B(k);
x(j)=x(j)+t(j);
end
x(j)=x(j)/k;
end
五、运行结果
[A,x]=Guass1(2)
A =
0.8889 0.5556 0.5556
x =
0 0.7746 -0.7746
>>
>> [A,x]=Guass1(1)
A =
1 1
xHale Waihona Puke =0.5774 -0.5774
西北农林科技大学数值分析数值法实验报告
数值法实验报告专业班级:信息与计算科学121 姓名:金辉 学号:20120142801)实验目的本次实验的目的是熟练《数值分析》第二章“插值法”的相关容,掌握三种插值法:牛顿多项式插值,三次样条插值,拉格朗日插值,并比较三种插值法的优劣。
本次试验要求编写牛顿多项式插值,三次样条插值,拉格朗日插值的程序编码,并在MATLAB 软件中去实现。
2)实验题目实验一:已知函数在下列各点的值为试用4次牛顿插值多项式P 4(x )及三次样条函数S (x )(自然边界条件)对数据进行插值。
用图给出{(x i ,y i ),x i =0.2+0.08i ,i=0,1, 11, 10},P 4(x )及S (x )。
实验二:在区间[-1,1]上分别取10,20n =用两组等距节点对龙格函数21()125f x x =+作多项式插值及三次样条插值,对每个n 值,分别画出插值函数即()f x 的图形。
实验三:可以得到平根函数的近似,在区间[0,64]上作图。
(1)用这9各点作8次多项式插值L8(x).(2)用三次样条(自然边界条件)程序求S(x)。
从结果看在[0,64]上,那个插值更精确;在区间[0,1]上,两种哪个更精确?3)实验原理与理论基础《数值分析》第二章“插值法”的相关容,包括:牛顿多项式插值,三次样条插值,拉格朗日4)实验容实验一:已知函数在下列各点的值为试用4次牛顿插值多项式P4(x)及三次样条函数S(x)(自然边界条件)对数据进行插值。
用图给出{(x i,y i),x i=0.2+0.08i,i=0,1, 11, 10},P4(x)及S(x)。
(1)首先我们先求牛顿插值多项式,此处要用4次牛顿插值多项式处理数据。
已知n次牛顿插值多项式如下:P n=f(x0)+f[x0,x1](x-x0)+ f[x0,x1,x2](x-x0) (x-x1)+···+f[x0,x1,···x n](x-x0) ···(x-x n-1)我们要知道牛顿插值多项式的系数,即均差表中得部分均差。
数值分析实验报告高斯消元法和列主消元法
《计算方法》实验指导书 实验三、高斯消元法和列主消元法一、实验目的:1. 通过matlab 编程解决高斯消元发和列主消元发来解方程组的问题, 加强编程能力和编程技巧,要熟练应用matlab 程序来解题,练习从数值分析的角度看问题进而来解决问题。
更深一步体会这门课的重要性,练习动手能力,同时要加深对数值问题的理解,要熟悉matlab 编程环境。
二、实验要求:用matlab 编写代码并运行高斯消元法和列主消元发来解下面的方程组的问题,并算出结果。
三、实验内容:用高斯消元法和列主消元法来解题。
1.实验题目:用高斯消元法和列主消元法来解下列线性方程组。
⎪⎪⎩⎪⎪⎨⎧−=+−−−=+−−=+−−=−+−.142,16422,0,13143214321432432x x x x x x x x x x x x x x x 2.实验原理高斯消元法:就是把方程组变成上三角型或下三角形的解法。
上三角形是从下往上求解,下三角形是从上向下求解,进而求得结果。
而列主消元法是和高斯消元法相类似,只不过是在开始的时候找出x1的系数的最大值放在方程组的第一行,再化三角形再求解。
3.设计思想高斯消元法:先把方程组的第一行保留,再利用第一行的方程将其余几行的含有x1的项都消去,再保留第二行,同理利用第二行的方程把第二行以下的几行的含有x2项的都消去,以此类推。
直到最后一行只含有一个未知数,化为上三角形,求得最后一行的这个未知数的值,再回带到倒数第二个方程求出另一个解,再依次往上回带即可求出这个方程组的值。
而列主消元法与高斯消元法类似,只不过在最开始时找出x1项系数的最大值与第一行交换再进行与高斯算法相似的运算来求出方程组的解。
4.源代码高斯消元法的程序:f unction [RA,RB,n,X]=gaus(A,b)B=[A b]; n=length(b); RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意:因为RA~=RB,所以此方程组无解.')returnendif RA==RBif RA==ndisp('请注意:因为RA=RB=n,所以此方程组有唯一X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1for k=p+1:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意:因为RA=RB<n,所以此方程组有无穷多解.')endend在工作窗口输入程序:A=[1 -1 1 -3; 0 -1 -1 1;2 -2 -4 6;1 -2 -4 1];b=[1;0; -1;-1]; [RA,RB,n,X] =gaus (A,b)请注意:因为RA=RB=n,所以此方程组有唯一解.运行结果为:RA =4RB =4n =4X =-0.50000.5000.列主消元发的程序:function [RA,RB,n,X]=liezhu(A,b)B=[A b]; n=length(b); RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意:因为RA~=RB,所以此方程组无解.')returnendif RA==RBif RA==ndisp('请注意:因为RA=RB=n,所以此方程组有唯一解.')X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1[Y,j]=max(abs(B(p:n,p))); C=B(p,:);B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;for k=p+1:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意:因为RA=RB<n,所以此方程组有无穷多解.')endend在工作窗口输入程序:A=[1 -1 1 -3; 0 -1 -1 1;2 -2 -4 6;1 -2 -4 1];b=[1;0; -1;-1]; [RA,RB,n,X]=liezhu(A,b)请注意:因为RA=RB=n,所以此方程组有唯一解.运行结果为:RA =4RB =4n =4X =-0.50000.5000实验体会:通过这次实验我了解了高斯消元法和列主消元方法的基本思想,虽然这两个程序的编写是有点困难的,但运行起来还是比较容易的,解决了不少实际问题的计算。
数值分析—实验报告1
wc=yi-jqz
wc = 0.00487093023460 数据分析:
表:ln(0.6)
内插
外插
x
0.4
0.5
0.7
0.7
0.8
0.9
f(x) -0.916291 -0.693147 -0.356675 -0.356675 -0.223144 -0.105361
yi
-0.50660833333333
-0.693147
-0.356675
试用 Lagrange 插值多项式来计算 ln(0.6)的近似值并估计误差。
运行结果:
format long
x=[0.40 0.50 0.70];y=[-0.916291 -0.693147
-0.356675];xin=0.6;yi=lg201541110131(x,y,xin)
实验仪器:
1、支持 Intel Pentium Ⅲ及其以上 CPU,内存 256MB 以上、硬盘 1GB 以上容量的微机; 软件配
有 Windows98/2000/XP 操作系统及 MATLAB 软件等。
2、了解 MATLAB 等软件的特点及系统组成,在电脑上操作 MATLAB 等软件。
实验内容、步骤及程序: 一、Lagrange 插值函数 程序: function yi=lg201541110131(x,y,xin) n=length(x); p=zeros(1,n); for k=1:n
一阶均差
2.23144000000000 -1.83026666666667
1.68236000000000
0
0
0
二阶均差
三阶均差
-0.916291000000000
数值分析实验——数值积分
桂林电子科技大学数学与计算科学学院实验报告 实验室: 06406 实验日期: 2014 年 11 月 21 日院(系) 数学与应用数学 年级、专业、班 姓名 成绩课程名称 数值分析实验 实验项目 名 称 实验积分 指导 教师李光云 一 、实验目的通过实验掌握利用Matlab 进行数值积分的操作,掌握Matlab 中的几种内置求积分函数,进一步理 解复化梯形,复化辛普生公式,并编程实现求数值积分二、实验原理Matlab 中,有内置函数计算积分:>> z = trapz(x,y)其中,输入x ,y 分别为已知数据的自变量和因变量构成的向量,输出为积分值。
>> z = quad(fun,a,b)这个命令是使用自适应求积的方法计算积分的命令。
其中,fun 为被积函数,a ,b 为积分区间。
我们还可以利用复化梯形公式()()()()⎪⎭⎫ ⎝⎛++-=∑-=11022n i i n x f x f x f n a b I 三、使用仪器,材料电脑 MATLAB四、实验内容与步骤1. 编写复化辛普生公式的Matlab 的程序。
2. 利用复化梯形法程序计算12041I dx x =+⎰,记录下计算结果随着n 增加的变化情况,画图与复化梯形公式的情况比较收敛速度。
3. 积分⎰dx x x sin 的原函数无法用初等函数表达,结合Matlab 复化梯形程序,用描点法绘制其原函数⎰x dt t t 1sin 在区间[]50,1的图形。
五、实验过程原始记录(数据,图表,计算等)一、复化Simpson公式程序:function s=Simpson(a,b,n)%输出s为积分的数值解,输入(a,b)为积分区间,n为等分区间的个数. h=(b-a)/(n*2);s1=0;s2=0;s=h*(f(a)+f(b))/3;%先计算特殊两点相加.for k=1:nx1=a+h*(2*k-1); %利用循环计算其他点的相加.s1=s1+f(x1);endfor k=1:(n-1)x2=a+h*2*k; %利用循环计算其他点的相加.s2=s2+f(x2);ends=s+h*(4*s1+2*s2)/3;画图程序format long;% k为等分区间个数,t存储积分值.k=2:1:40;for i=1:length(k)t(i)=Simpson(0,1,k(i));disp([k(i),t(i)]);endplot(k,t,'.','MarkerSize',20)二、Simpson(0,1,26)ans =3.14159265358779>> Untitled62.000000000000003.141568627450983.00000000000000 3.141591780936044.00000000000000 3.141592502458715.00000000000000 3.141592613939226.00000000000000 3.141592640305387.00000000000000 3.141592648320658.00000000000000 3.141592651224829.00000000000000 3.1415926524231710.00000000000000 3.1415926529697911.00000000000000 3.1415926532398112.00000000000000 3.1415926533821513.00000000000000 3.1415926534613414.00000000000000 3.1415926535074515.00000000000000 3.1415926535353616.00000000000000 3.1415926535528417.00000000000000 3.1415926535641118.00000000000000 3.1415926535715619.00000000000000 3.1415926535766120.00000000000000 3.1415926535801121.00000000000000 3.1415926535825622.00000000000000 3.1415926535843223.00000000000000 3.1415926535856024.00000000000000 3.1415926535865525.00000000000000 3.1415926535872526.00000000000000 3.1415926535877927.00000000000000 3.1415926535881928.00000000000000 3.1415926535885129.00000000000000 3.1415926535887530.00000000000000 3.1415926535889431.00000000000000 3.1415926535890932.00000000000000 3.1415926535892233.00000000000000 3.1415926535893134.00000000000000 3.1415926535893935.00000000000000 3.1415926535894636.00000000000000 3.1415926535895137.00000000000000 3.1415926535895538.00000000000000 3.1415926535895939.00000000000000 3.1415926535896240.00000000000000 3.14159265358964由图形可以知道,复化Simpson求积公式的收敛速度比复化梯形求积公式收敛速度快三、b=2:1:50;>> y=tixing(1,b,26)Warning: Colon operands must be real scalars.> In tixing at 3y =Columns 1 through 40.65931335927431 1.31862671854861 1.97794007782292 2.63725343709723 Columns 5 through 83.29656679637153 3.955880155645844.615193514920155.27450687419446 Columns 9 through 125.933820233468766.593133592743077.25244695201738 7.91176031129169 Columns 13 through 168.57107367056599 9.23038702984030 9.88970038911461 10.54901374838891 Columns 17 through 2011.20832710766322 11.86764046693753 12.52695382621183 13.18626718548614 Columns 21 through 2413.84558054476045 14.50489390403476 15.16420726330906 15.82352062258337 Columns 25 through 2816.48283398185768 17.14214734113198 17.80146070040629 18.46077405968060 Columns 29 through 3219.12008741895491 19.77940077822921 20.43871413750352 21.09802749677782 Columns 33 through 3621.75734085605213 22.41665421532644 23.07596757460075 23.73528093387506 Columns 37 through 4024.39459429314936 25.05390765242366 25.71322101169798 26.37253437097228 Columns 41 through 4427.03184773024659 27.69116108952090 28.35047444879520 29.00978780806952Columns 45 through 4829.66910116734383 30.32841452661813 30.98772788589244 31.64704124516674Column 4932.30635460444105>> x=1:1:49;>> plot(x,y)六、实验结果分析或总结通过这次实验,学会了如何使用Matlab进行数值积分的操作以及几种内置求积分函数的方法,但编程一直都是自己的弱势,一些程序的代码还是无法自己编译出来。
