研究生数值分析上机试题及解答
- 格式:doc
- 大小:60.00 KB
- 文档页数:4
东华大学研究生数值分析试题(上机部分)A 卷2008年12月 时间:60分钟班级 学号 机号 姓名 得分 注意:要求写出M 函数(如果需要)、MATLAB 命令和计算结果。
1. 求下列方程组在0<α, β<1中的解⎩⎨⎧-=+=βαββααsin 2.0cos 7.0cos 2.0sin 7.0 命令fun=inline('[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))]','x'); [x,f,h]=fsolve(fun,[0.5 0.5]) 结果α=0.5265,β=0.50792命令>> fun=inline('c(1)+c(2)*x.^2','c','x'); >> x=[1.1 1.3 1.4 1.6 1.8]; >> y=[26 22 23 24 25];>> c=lsqcurvefit(fun,[0 0],x,y) 结果 c =23.7256 0.12873.求解下列微分方程组2(0)2013(0)1x x yx t y x yy '=-=⎧<<⎨'=+=⎩(结果只要求写出t =1时的解) 命令>> fun=inline('[y(1)-2*y(2);3*y(1)+y(2)]','t','y'); >> [t,y]=ode45(fun,[0 1], [2 1]) 结果x(1)=-5.6020, y(1)=2.15634.用定步长Gauss 积分法(课本123页)计算积分31e ln(1)x x dx -+⎰的近似值(等分数取4,每段取2个Gauss 点)。
命令fun=inline('exp(-x).*log(1+x)','x'); nagsint(fun,1,3,4,2) 结果 0.30865.矩阵改进平方根分解(课本25页)的计算公式为: d 1=a 11, 对i =2, 3, ⋯, n ,iki k ik ii i j ij ij j k jk ik ij ij l s a d i j d s l l s a s ∑∑-=-=-=-==-=1111,1,,2,1 ,/ ,试编写矩阵改进平方根分解的程序,并求矩阵1111551514A -⎛⎫ ⎪=-- ⎪ ⎪-⎝⎭的改进平方根分解。
数值分析上机题目4(总21页) --本页仅作为文档封面,使用时请直接删除即可----内页可以根据需求调整合适字体及大小--实验一实验项目:共轭梯度法求解对称正定的线性方程组 实验内容:用共轭梯度法求解下面方程组(1) 123421003131020141100155x x x x -⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪=⎪ ⎪ ⎪-- ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ 迭代20次或满足()(1)1110k k x x --∞-<时停止计算。
编制程序:储存m 文件function [x,k]=CGmethod(A,b)n=length(A);x=2*ones(n,1);r=b-A*x;rho=r'*r; k=0;while rho>10^(-11) & k<1000 k=k+1; if k==1 p=r; elsebeta=rho/rho1; p=r+beta*p; end w=A*p;alpha=rho/(p'*w); x=x+alpha*p; r=r-alpha*w; rho1=rho;rho=r'*r; end运行程序: clear,clcA=[2 -1 0 0;-1 3 -1 0;0 -1 4 -1;0 0 -1 5]; b=[3 -2 1 5]'; [x,k]=CGmethod(A,b)运行结果: x =(2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵, A[i,i]=4,A[i,i-1]=A[i-1,i]=-1,i=2,3,..,n b[1]=3, b[n]=3, b[i]=2,i=2,3,…,n-1迭代10000次或满足()()710k k r b Ax -=-≤时停止计算。
编制程序:储存m 文件function [x,k]=CGmethod_1(A,b) n=length(A);x(1:n,1)=0;r=b-A*x;r1=r; k=0;while norm(r1,1)>10^(-7)&k<10^4 k=k+1; if k==1 p=r; elsebeta=(r1'*r1)/(r'*r);p=r1+beta*p; end r=r1; w=A*p;alpha=(r'*r)/(p'*w); x=x+alpha*p; r1=r-alpha*w; end运行程序: clear,clc n=1000; A=hilb(n); b=sum(A')';[x,k]=CGmethod(A,b)实验二1、 实验目的:用复化Simpson 方法、自适应复化梯形方法和Romberg 方法求数值积分。
源程序如下:function[x_jac,x_gauss,x_sor,A,b]=diedai(n,w,e,N,Ain,bin) if nargin<3e=1e-5;N=1000;endif nargin<4N=1000;endif nargin<2error('²ÎÊý²»¹»');endA=hilb(n);b=A*ones(n,1);if nargin>4if nargin==5error('ÇëÊäÈëb');elseA=Ain;b=bin;endend%%Ñſ˱ȵü´údiedai_n=0;x=zeros(n,2);while max(abs(b-A*x(:,2)))>ediedai_n=diedai_n+1;if diedai_n>Ndisp('Ñſ˱ȵü´úδ´ïµ½¾«¶È');break;endx(:,1)=x(:,2);for i=1:nif i>1x(i,2)=(b(i)-A(i,1:(i-1))*x(1:(i-1),1)-A(i,(i+1):n)*x((i+ 1):n,1))/A(i,i);elsex(i,2)=(b(i)-A(i,(i+1):n)*x((i+1):n,1))/A(i,i);endendendx_jac.x=x(:,2);x_jac.n=diedai_n-1;x_jac.error=max(abs(b-A*x(:,2)));%%¸ß˹µü´údiedai_n=0;x=zeros(n,2);while max(abs(b-A*x(:,2)))>ex(:,1)=x(:,2);diedai_n=diedai_n+1;if diedai_n>Ndisp('¸ß˹¡ªÈüµÂ¶ûµü´úδ´ïµ½¾«¶È');break;endfor i=1:nif i>1x(i,2)=(b(i)-A(i,1:(i-1))*x(1:(i-1),2)-A(i,(i+1):n)*x((i+ 1):n,1))/A(i,i);elsex(i,2)=(b(i)-A(i,(i+1):n)*x((i+1):n,1))/A(i,i); endendendx_gauss.x=x(:,2);x_gauss.n=diedai_n-1;x_gauss.error=max(abs(b-A*x(:,2)));%%SORµü´údiedai_n=0;x=zeros(n,2);while max(abs(b-A*x(:,2)))>ex(:,1)=x(:,2);diedai_n=diedai_n+1;if diedai_n>Ndisp('SORµü´úδ´ïµ½¾«¶È');break;endfor i=1:nif i>1x(i,2)=x(i,1)+w*(b(i)-A(i,1:(i-1))*x(1:(i-1),2)-A(i,i:n)* x(i:n,1))/A(i,i);elsex(i,2)=x(i,1)+w*(b(i)-A(i,i:n)*x(i:n,1))/A(i,i);endendendx_sor.x=x(:,2);x_sor.n=diedai_n-1;x_sor.error=max(abs(b-A*x(:,2)));end运行结果:由于SOR的松弛因子取1时与高斯-赛德尔迭代相同,故在矩阵n取8,10时省去了SOR迭代松弛因子取1的步骤。
------------------------------------------------ 装 ---------------------------------订 ---------------------------------线 ------------------------------------------------装 订 线 左 侧 不 要 书 写 内 容允许使用计算器一、 填空题 (本大题共10小题,每小题 2分,共 20分)1. 若2.71828x e == ,取近似值* 2.7180x =,则*x 具有 4 位有效数字。
2.为了提高数值计算精度,应将8格式进行计算。
3.已知n=3时牛顿—柯特斯系数(3)(3)(3)012133,,888C C C ===,那么(3)3C =18 。
4.设3()1f x x x =+-,则函数的四阶差商[0,1,2,3,4]f = 0 。
5. 用牛顿迭代法解方程0x x e --=在0.5x =附近的近似实根的牛顿迭代格式为)1,0(e 1e )()(1=+--='-=--+n x x x f x f x x nnx x n n n n n n6. 对给定的剖分01:n a x x x b ∆=<<<= ,当()s x 满足条件 ()s x 在[a,b]有2阶连续导数且在每个子区间上是个3次多项式 时是三次样条函数。
7.用最小二乘法拟合三点()()()0,1,1,3,2,2A B C 的直线是1322y x =+。
8.向量序列()211cos ,sin ,3Tk k x e k k k k -⎛⎫=+ ⎪⎝⎭ 的极限向量为()0,1,3T9.求积公式 10311()()(1)434f x dx f f ≈+⎰的代数精度为 2 。
10.若绝对误差限为31102-⨯,那么近似数0.03600有 2 位有效数字二、单项选择题(本大题共5小题,每小题 2 分,共 10分)1. 已知实验数据555521111(,)(1,2,3,4,5),15,31,55,105.5,k k k k kk k k k k k x y k x y x x y =========∑∑∑∑其中则用最小二乘法求近似公式01y a a x =+的法方程为( C )A 0101153155105.5a a a a +=⎧⎨+=⎩B 0101515551531105.5a a a a +=⎧⎨+=⎩C 0101515311555105.5a a a a +=⎧⎨+=⎩ D0101531153155105.5a a a a +=⎧⎨+=⎩ 2. 以下矩阵是严格对角占优矩阵的是( B )A 3210141011410012⎛⎫ ⎪ ⎪ ⎪⎪⎝⎭ B 2100131013610113-⎛⎫⎪--⎪ ⎪-- ⎪-⎝⎭C 5210113121410012-⎛⎫⎪--⎪ ⎪⎪⎝⎭D 4211141021411315⎛⎫⎪ ⎪⎪- ⎪⎝⎭3.已知两种递推公式11(1)35(1,2,,20)31(2)(20,,1)55n n n n I nI n I I n n n--=-==-= 则在数值计算过程中( C )。
数值分析上机题第一章17.(上机题)舍入误差与有效数 设∑=-=Nj N j S 2211,其精确值为)111-23(21+-N N 。
(1)编制按从大到小的顺序1-1···1-311-21222N S N +++=,计算N S 的通用程序;(2)编制按从小到大的顺序121···1)1(111222-++--+-=N N S N ,计算NS 的通用程序;(3)按两种顺序分别计算210S ,410S ,610S ,并指出有效位数(编制程序时用单精度); (4)通过本上机题,你明白了什么?解: 程序:(1)从大到小的顺序计算1-1···1-311-21222N S N +++=:function sn1=fromlarge(n) %从大到小计算sn1format long ; sn1=single(0); for m=2:1:nsn1=sn1+1/(m^2-1); end end(2)从小到大计算121···1)1(111222-++--+-=N N S N function sn2=fromsmall(n) %从小到大计算sn2format long ; sn2=single(0); for m=n:-1:2sn2=sn2+1/(m^2-1); end end(3)总的编程程序为: function p203()clear allformat long;n=input('please enter a number as the n:') sn=1/2*(3/2-1/n-1/(n+1));%精确值为snfprintf('精确值为%f\n',sn);sn1=fromlarge(n);fprintf('从大到小计算的值为%f\n',sn1);sn2=fromsmall(n);fprintf('从小到大计算的值为%f\n',sn2);function sn1=fromlarge(n) %从大到小计算sn1 format long;sn1=single(0);for m=2:1:nsn1=sn1+1/(m^2-1);endendfunction sn2=fromsmall(n) %从小到大计算sn2 format long;sn2=single(0);for m=n:-1:2sn2=sn2+1/(m^2-1);endendend运行结果:从而可以得到N值真值顺序值有效位数2 100.740050 从大到小0.740049 5从小到大0.740050 64 100.749900 从大到小0.749852 3从小到大0.749900 66 100.749999 从大到小0.749852 3从小到大0.749999 6(4)感想:通过本上机题,我明白了,从小到大计算数值的精确位数比较高而且与真值较为接近,而从大到小计算数值的精确位数比较低。
数值分析考试及答案作者:日期:班级• • •• • •• • •• • • o • • •学号• • •• • •姓名密• • •• • •o• • •• • •东北大学研究生院考试试卷2011 —2012 学年第一学期课程名称:数值分析(共3页)一、解答下列各题:(每题5分,共30分)1.设近似值x具有5位有效数字,则x的相对误差限为多少? 解:记x* 0.吋2…10m,则x的相对误差为:0.5 10m 50.a1a2... 10m0.5 10 50.10.5 10即,相对误差限为:0.5 102.问a, b满足什么条件时,矩阵Ao • • •• • •• • •线总分一——二三四五4 2 02 5a有分解式A GG T,并求a b 2时0 b 54 2 0 2 1 0解:由于A 2 5 a 1 2 a/2 (A对称正定时)0 b 5 0 b/2 5 ab/4所以,当2 .5 a b 2 5时有分解式 A GG T,a b 2 时有:4 2 0 2 0 0 2 1 0A 2 5 2 1 2 0 0 2 10 2 5 0 1 2 0 0 23.解线性方程组X1 2x2 2 的Jacobi 迭代法是否收敛,为什么?2x19x2 3的分解式(其中G是对角线元素大于零的下三角形矩阵)解:Jacobi迭代矩阵为:B2/92,所以,(B) 2/3 1所以,Jacobi迭代法是否收敛.4.对方程f (x) (x3 a)20建立敛?若收敛,收敛阶是多少?解:Newton迭代格式为:X k 1 X kf(xk)f (X k)由于迭代函数为:(x)?X ka6x2所以,此迭代格式收敛,收敛阶是Newton迭代格式,并说明此迭代格式是否收3X k a2~ ,x k6X k,方程根为:1.56k 6:k2, k 012-3 a,所以,5.设f (x) 4x3 3x 5,求差商f[0,1], f[1,2,3,4]和f [1,2,3,4,5]。
● 当22()3x x ϕ+=时,'12()3x x ϕ=,因此'1(2) 1.3333ϕ=>1,,因此,该迭代格式不收敛。
● 当2()x ϕ='2()x ϕ=,因此'2(2)0.75ϕ=<1,,因此,该迭代格式收敛。
● 当32()3x x ϕ=-时,'322()x xϕ=,因此'3(2)0.5ϕ=<1,,因此,该迭代格式收敛。
● 当242()23x x x ϕ-=- 时,2'44224()1213x x x x x ϕ--=-+,因此'4(2)0ϕ=<1,,因此,该迭代格式收敛。
(2)、● 当22()3x x ϕ+=时,迭代法计算公式是20122.5,3k k x x x ++==,程序如下: >> fi=inline('(x.*x+2)/3');x0=2.5;er=1;k=0;while er>0.00001x=fi(x0);er=abs(x-x0);x0=xk=k+1end运行结果如下:x0 =2.7500k =1x0 =3.1875k =2x0 =4.0534k =3x0 =6.1433k =x0 =13.2468k =5x0 =59.1589k =6x0 =1.1673e+003 k =7x0 =4.5416e+005 k =8x0 =6.8755e+010 k =9x0 =1.5757e+021 k =10x0 =8.2765e+041 k =11x0 =2.2834e+083 k =12x0 =1.7379e+166 k =13x0 =Infk =14x0 =Infk =15由以上计算结果看,序列是发散的,运行14次已经超出计算机的识别范围,当2()x ϕ迭代法计算公式是1k x +=程序运行结果如下:>> fi=inline('sqrt(3*x-2)');x0=2.5;er=1;k=0;while er>0.00001x=fi(x0);er=abs(x-x0);x0=xk=k+1;endx0 =2.3452x0 =2.2440x0 =2.1753x0 =2.1274x0 =2.0934x0 =2.0689x0 =2.0510x0 =2.0379x0 =2.0282x0 =2.0211x0 =2.0157x0 =2.0118x0 =2.0088x0 =2.0066x0 =2.0049x0 =2.0037x0 =2.0028x0 =2.0021x0 =2.0016x0 =2.0012x0 =2.0009x0 =2.0007x0 =2.0005x0 =2.0004x0 =2.0003x0 =2.0002x0 =2.0002x0 =2.0001x0 =2.0001x0 =2.0001x0 =2.0000x0 =2.0000x0 =2.0000>>由以上计算结果看,序列收敛与2,所以x=2是f(x)= 232x x -+=0的根。
如有帮助欢迎下载支持数值分析上机题姓名:陈作添学号: 040816习题 120.(上机题)舍入误差与有效数N11 3 1 1设S N,其精确值为 。
22 2 NN 1j 2j1(1)编制按从大到小的顺序111 ,计算 S 的通用程序。
S N1 321N 21 N22(2)编制按从小到大的顺序111,计算 S 的通用程序。
S N1(N 1)2 122 1NN 2(3)按两种顺序分别计算S 102 , S 104 , S 106 ,并指出有效位数。
(编制程序时用单精度)(4)通过本上机题,你明白了什么?按从大到小的顺序计算 S N 的通用程序为: 按从小到大的顺序计算 S N 的通用程序为:#include<iostream.h> #include<iostream.h> float sum(float N) float sum(float N) {{float j,s,sum=0; float j,s,sum=0; for(j=2;j<=N;j++) for(j=N;j>=2;j--) {{s=1/(j*j-1); s=1/(j*j-1); sum+=s;sum+=s;}}return sum;return sum;}}从大到小的顺序的值从小到大的顺序的值精确值有效位数从大到小从小到大0.7400490.740050.74004965 S 1020.7498520.74990.74994 4S 1040.7498520.7499990.74999936S 106通过本上机题, 看出按两种不同的顺序计算的结果是不相同的,按从大到小的顺序计算的值与精确值有较大的误差, 而按从小到大的顺序计算的值与精确值吻合。
从大到小的顺序计算得到的结果的有效位数少。
计算机在进行数值计算时会出现“大数吃小数”的现象,导致计算结果的精度有所降低,我们在计算机中进行同号数的加法时,采用绝对值较小者先加的算法,其结果的相对误差较小。
++中的待定系数,使其A f(1)(0)武汉理工大学研究生课程考试标准答案用纸课程名称:数值计算(A ) 任课教师 :一. 简答题,请简要写出答题过程(每小题5分,共30分) 3.14159265358979的近似值,它们各有几位有效数字,绝对误差和相对误差分别是多少?3分)2分)2.已知()8532f x x x =+-,求0183,3,,3f ⎡⎤⎣⎦,0193,3,,3f ⎡⎤⎣⎦.(5分)3.确定求积公式10120()(0)(1)(0)f x dx A f A f A f '≈++⎰中的待定系数,使其代数精度尽量高,并指明该求积公式所具有的代数精度。
解:要使其代数精度尽可能的高,只需令()1,,,m f x x x =使积分公式对尽可能大的正整数m 准确成立。
由于有三个待定系数,可以满足三个方程,即2m =。
由()1f x =数值积分准确成立得:011A A += 由()f x x =数值积分准确成立得:121/2A A += 由2()f x x =数值积分准确成立得:11/3A =解得1201/3,1/6,2/3.A A A === (3分)此时,取3()f x x =积分准确值为1/4,而数值积分为11/31/4,A =≠所以该求积公式的最高代数精度为2次。
(2分)4.求矩阵101010202A -⎡⎤⎢⎥=⎢⎥⎢⎥-⎣⎦的谱半径。
解 ()()101101322I A λλλλλλλ--=-=--- 矩阵A 的特征值为1230,1,3λλλ=== 所以谱半径(){}max 0,1,33A ρ== (5分)5. 设10099,9998A ⎛⎫= ⎪⎝⎭计算A 的条件数()(),2,p cond A P =∞.解:**19899-98999910099-100A A A A --⎛⎫⎛⎫=⇒== ⎪ ⎪-⎝⎭⎝⎭矩阵A 的较大特征值为198.00505035,较小的特征值为-0.00505035,则1222()198.00505035/0.0050503539206cond A A A -=⨯==(2分)1()199********c o n d A A A -∞∞∞=⨯=⨯=(3分)22001130101011010220100110110()(12)()(12)()()()()()x x x x x x x x H x y y x x x x x x x x x x x x x x y x x y x x x x ----=-+-------''+-+---(5分)并依条件1(0)1,(0),(1)2,(1) 2.2H H H H ''====,得2222331()(12)(1)2(32)(1)2(1)211122H x x x x x x x x x x x =+-+-+-+-=++ (5分)2.已知()()()12,11,21f f f -===,求()f x 的Lagrange 插值多项式。
东华大学研究生数值分析试题(上机部分)
A 卷2008年12月 时间:60分钟
班级 学号 机号 姓名 得分 注意:要求写出M 函数(如果需要)、MATLAB 命令和计算结果。
1. 求下列方程组在0<α, β<1中的解
⎩
⎨
⎧-=+=βαββ
ααsin 2.0cos 7.0cos 2.0sin 7.0 命令
fun=inline('[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))]','x'); [x,f,h]=fsolve(fun,[0.5 0.5]) 结果
α=0.5265,β=0.5079
2
命令
>> fun=inline('c(1)+c(2)*x.^2','c','x'); >> x=[1.1 1.3 1.4 1.6 1.8]; >> y=[26 22 23 24 25];
>> c=lsqcurvefit(fun,[0 0],x,y) 结果 c =
23.7256 0.1287
3.求解下列微分方程组2(0)2013(0)1
x x y
x t y x y
y '=-=⎧<<⎨
'=+=⎩
(结果只要求写出t =1时的解) 命令
>> fun=inline('[y(1)-2*y(2);3*y(1)+y(2)]','t','y'); >> [t,y]=ode45(fun,[0 1], [2 1]) 结果
x(1)=-5.6020, y(1)=2.1563
4.用定步长Gauss 积分法(课本123页)计算积分3
1e ln(1)x x dx -+⎰的近似值(等
分数取4,每段取2个Gauss 点)。
命令
fun=inline('exp(-x).*log(1+x)','x'); nagsint(fun,1,3,4,2) 结果 0.3086
5.矩阵改进平方根分解(课本25页)的计算公式为: d 1=a 11, 对i =2, 3, ⋯, n ,
ik
i k ik ii i j ij ij j k jk ik ij ij l s a d i j d s l l s a s ∑∑-=-=-=-==-=1
1
1
1,
1,,2,1 ,/ ,
试编写矩阵改进平方根分解的程序,并求矩阵1111551514A -⎛⎫ ⎪
=-- ⎪ ⎪-⎝⎭
的改进
平方根分解。
%M 函数
function [l,d]=ldlt(a) n=length(a);
l=zeros(n,n);s=l;d=zeros(1,n); d(1)=a(1,1); for i=2:n
for j=1:i-1
t=0;for k=1:j-1,t=t+s(i,k)*l(j,k);end; s(i,j)=a(i,j)-t;l(i,j)=s(i,j)/d(j); end;
t=0;for k=1:i-1,t=t+s(i,k)*l(i,k);end;d(i)=a(i,i)-t; end
l=eye(n)+l; 命令
>> [l,d]=ldlt([1 -1 1;-1 5 -5;1 -5 14]) 结果
l = 1 0 0 -1 1 0 1 -1 1 d = 1 4 9
东华大学研究生数值分析试题(上机部分)
B 卷2008年12月 时间:60分钟
班级 学号 机号 姓名 得分 注意:要求写出M 函数(如果需要)、MATLAB 命令和计算结果。
1.求积分1
0sin(2)
x
x dx e
⎰的近似值 命令
>> fun=inline('sin(2*x)./exp(x)','x'); >> quadl(fun, 0,1) 结果 0.3943
2.求下列方程组在0<x , y , z <2中的解
⎪⎩
⎪⎨⎧=+=-+=--8101110741222
2z y z y x z y x 命令 >>
fun=inline('[12*x(1)-x(2)^2-4*x(3)-7;x(1)^2+10*x(2)-x(3)-11;x(2)^2+10*x(3)-8]'); >> [x,f,h]=fsolve(fun,[1 1 1]) 结果
x = 0.9089, y=1.0856, z=0.6821 3.已知函数。
命令
>> x=[1 2 3 4];y=[4.6 4 3.3 3.4]; >> pp=csape(x,y,'complete',[0 0]); >> pp.coefs 结果
s(x)= 0.28(x-1)3 -0.88(x-1)2+4.6, 1≤x ≤2 0.26(x-1)3 - 0.04(x-1)2 -0.92(x-1)+4, 2≤x ≤3 -0.42(x-1)3 +0.74(x-1)2-0.22 (x-1)+3.3, 3≤x ≤4
4.用课本157页改进Euler 法程序(h=0.2)解下列微分方程
32 1
)2(ln '<<⎩⎨
⎧=+=x y y
x y
(结果只要求写出x =3时y 的解) 命令
>> fun=inline('log(y(1))+y(2)','x','y'); >> [t,y]=naeuler2s(fun,[2 3],1,0.2) 结果
y(3)=4.1899
5.矩阵Crout 分解(课本23页)的计算公式为,对i =1, 2, ⋯, n ,
n
i i j l u l a u n i i j u l a l ii i k kj ik ij ij i k ki jk ji ji ,,2,1,/)(,
,,1,,1
1
1
1
++=-=+=-=∑∑-=-=
试编写矩阵Crout 分解的程序,并求矩阵120104107A ⎛⎫ ⎪
=-- ⎪ ⎪⎝⎭
的Crout 分解。
%M 函数
function [l,u]=crout(a)
n=length(a); l=zeros(n,n);u=l; for i=1:n
j=i;t=0;for k=1:i-1,t=t+l(j,k)*u(k,i);end;l(j,i)=a(j,i)-t; for j=i+1:n
t=0;for k=1:i-1,t=t+l(j,k)*u(k,i);end;l(j,i)=a(j,i)-t;
t=0;for k=1:i-1,t=t+l(i,k)*u(k,j);end;u(i,j)=(a(i,j)-t)/l(i,i); end end
u=eye(n)+u;
命令
>> [l,u]=crout([1 2 0;-1 0 -4;1 0 7])
结果
l = 1 0 0 -1 2 0 1 -2 3 u = 1 2 0 0 1 -2 0 0 1。