大连理工大学-矩阵大作业

（1）从大到小的顺序的计算程序：

function y=snd(n) format long y=0; if n

disp('请输入大于1的数！') else s=0; i=2;

while i

s=single(s+(1/(i^2-1))); i=i+1; end y=s; end

（2）从小到大的顺序的计算程序：

function y=snx(n) format long y=0; if n

disp('请输入大于1的数！') else s=0; i=n; while 1

s=single(s+(1/(i^2-1))); i=i-1; if i==1 break end end

y=s; end

（3）按两种顺序分别计算① S 102的计算结果：

② S 104的计算结果：

③ S 106的计算结果：

计算时的有效位数为七位数。

并指出有效位数（编制程序时用单精度）

① 秦九昭算法计算程序：

function y=qjz(a,x) j=3;

i=size(a,2); switch i case 1 y=a(1); case 2

y=a(1)*x+a(2); otherwise

p=a(1)*x+a(2); while j

② 计算

计算结果如下：

在点23的值。

当x =23时f (x )=86652。

① Gauss 法计算程序和结果：程序：

A(1,:)=[31,-13,0,0,0,-10,0,0,0]; A(2,:)=[-13,35,-9,0,-11,0,0,0,0]; A(3,:)=[0,-9,31,-10,0,0,0,0,0]; A(4,:)=[0,0,-10,79,-30,0,0,0,-9]; A(5,:)=[0,0,0,-30,57,-7,0,-5,0]; A(6,:)=[0,0,0,0,-7,47,-30,0,0]; A(7,:)=[0,0,0,0,0,-30,41,0,0]; A(8,:)=[0,0,0,0,-5,0,0,27,-2]; A(9,:)=[0,0,0,-9,0,0,0,-2,29]; B=[-15;27;-23;0;-20;12;-7;7;10]; [a,b]=gauss(A,B); j=size(a,2); while j>=1 k=j+1; s=b(j); while k

s=s-x(k)*a(j,k); k=k+1; end

x(j)=s/a(j,j); j=j-1; end disp(x)

function [x,y]=gauss(a,b) num_i=size(a,1); j=1;

while j

while i

b(i,:)=b(i,:)-r*b(j,:); i=i+1; end j=j+1; end x=a; y=b;

运行的结果为：

x =(-0. 2890. 345-0. 713-0. 221-0. 4300. 154-0. 0580. 2010. 290)。

② 列主元消去法计算程序和结果：

A(1,:)=[31,-13,0,0,0,-10,0,0,0]; A(2,:)=[-13,35,-9,0,-11,0,0,0,0]; A(3,:)=[0,-9,31,-10,0,0,0,0,0]; A(4,:)=[0,0,-10,79,-30,0,0,0,-9]; A(5,:)=[0,0,0,-30,57,-7,0,-5,0]; A(6,:)=[0,0,0,0,-7,47,-30,0,0]; A(7,:)=[0,0,0,0,0,-30,41,0,0]; A(8,:)=[0,0,0,0,-5,0,0,27,-2]; A(9,:)=[0,0,0,-9,0,0,0,-2,29]; B=[-15;27;-23;0;-20;12;-7;7;10]; [a,b]=lzy(A,B); j=size(a,2); while j>=1 k=j+1; s=b(j); while k

s=s-x(k)*a(j,k); k=k+1; end

x(j)=s/a(j,j); j=j-1; end disp(x)

function [a1,b1]=lzy(a,b) [num_i,num_j]=size(a); ab=zeros(num_i,num_j+1); for k=1:num_j ab(:,k)=a(:,k); end

ab(:,num_j+1)=b(:,1); j=1;

while j

[max,max_i]=searmax(j,j,ab); i_nu=ab(j,:); ab(j,:)=ab(max_i,:); ab(max_i,:)=i_nu; m=j+1;

while m

ab(m,n)=ab(m,n)-(ab(m,j)/max)*ab(j,n); end m=m+1; end j=j+1; end

a1=zeros(num_i,num_j); for k=1:num_i a1(:,k)=ab(:,k); end

b1=ab(:,num_i+1);

function [b,c]=searmax(i,j,a) num_i=size(a,1); k=i;

m=abs(a(k,j)); c=k;

while k

if m>=abs(a(k,j)) continue else

m=abs(a(k,j)); c=k; end end b=a(c,j);

运行的结果为：

x =

(-0. 2370. 477-0. 767-0. 078-0. 3080. 146-0. 1710. 2850. 345)

（1）LU 分解的计算程序及结果：

function [l,u]=lufz(a) [num_i,num_j]=size(a);

l=eye(num_i); u=eye(num_i); for j=1:num_j u(1,j)=a(1,j); l(j,1)=a(j,1)/u(1,1); end i=2;

while i

while j

s=s+l(i,k)*u(k,j); end

u(i,j)=a(i,j)-s; s=0; for k=1:i-1

s=s+l(j+1,k)*u(k,i); end

l(j+1,i)=(a(j+1,i)-s)/u(i,i); j=j+1; end s=0; for k=1:i-1

s=s+l(i,k)*u(k,num_i); end

u(i,num_i)=a(i,num_i)-s; i=i+1; end

输入要求解的矩阵后运行的结果为：

⎛ 1 -0. 4191 0-0. 3051 00-0. 3541L = 000-0. 3981 0000-0. 1571 00000-0. 6551 0000-0. 112-0. 018-0. 0251 ⎝000-0. 119-0. 083-0. 014-0. 020-0. 092⎫⎪⎪⎪⎪

⎪

⎪⎪

⎪⎪⎪⎪⎭

29. 548-90-11 28. 259-10-3. 350

75. 461-31. 186

U = 44. 602

⎝

⎪⎪

-1. 277000⎪

⎪

-0. 45200-9⎪-7. 1800-5-3. 578⎪⎪45. 8732-30-0. 785-0. 562⎪

⎪

21. 381-0. 513-0. 367⎪

26. 413-2. 420⎪

⎪27. 390⎪⎭

-4. 194000

（2）带列主元的LU 分解计算程序及结果

由于Matlab 中自带带列主元的LU 分解函数，故计算程序如下：

运行结果如下：

⎛1⎫ ⎪

1 -0. 419⎪

0⎪-0. 3051 ⎪

0-0. 3541 0⎪

⎪

L = 000-0. 3981 ⎪

0⎪000-0. 1571 ⎪

00000-0. 6551 ⎪ 0⎪000-0. 112-0. 018-0. 0251 ⎪ 0⎪00-0. 119-0. 083-0. 014-0. 020-0. 0921⎝⎭

29. 548-90-11 28. 259-10-3. 350

75. 461-31. 186

U = 44. 602

⎪⎪

-1. 277000⎪

⎪

-0. 45200-9⎪-7. 1800-5-3. 578⎪⎪45. 8732-30-0. 785-0. 562⎪

⎪

21. 381-0. 513-0. 367⎪

-4. 194000

⎝⎛ 100000000⎫

010000000⎪

⎪ 001000000⎪

000100000⎪

⎪P = 0

00010000⎪ 000001000⎪为单位阵。 ⎪

000000100⎪

⎪ 000000010⎪

⎝

000000001⎪

⎪⎭

（3）求A 的逆矩阵的计算程序及结果：

A(1,:)=[31,-13,0,0,0,-10,0,0,0]; A(2,:)=[-13,35,-9,0,-11,0,0,0,0]; A(3,:)=[0,-9,31,-10,0,0,0,0,0]; A(4,:)=[0,0,-10,79,-30,0,0,0,-9]; A(5,:)=[0,0,0,-30,57,-7,0,-5,0]; A(6,:)=[0,0,0,0,-7,47,-30,0,0]; A(7,:)=[0,0,0,0,0,-30,41,0,0]; A(8,:)=[0,0,0,0,-5,0,0,27,-2]; A(9,:)=[0,0,0,-9,0,0,0,-2,29]; [num_i,num_j]=size(A); A_n=eye(num_i); E=eye(num_i); [L,U]=lufz(A); for i=1:num_i y=hd(1,L,E(:,i)); A_n(:,i)=hd(0,U,y); end disp(A_n)

function [l,u]=lufz(a) [num_i,num_j]=size(a);

26. 413-2. 420⎪

27. 390⎪⎪⎭

l=eye(num_i); u=eye(num_i); for j=1:num_j u(1,j)=a(1,j); l(j,1)=a(j,1)/u(1,1); end i=2;

while i

while j

s=s+l(i,k)*u(k,j); end

u(i,j)=a(i,j)-s; s=0; for k=1:i-1

s=s+l(j+1,k)*u(k,i); end

l(j+1,i)=(a(j+1,i)-s)/u(i,i); j=j+1; end s=0; for k=1:i-1

s=s+l(i,k)*u(k,num_i); end

u(i,num_i)=a(i,num_i)-s; i=i+1; end

function x=hd(f,a,b) [num_i,num_j]=size(a); x=zeros(num_i,1); switch f case 0 i=num_i-1;

x(num_i)=b(num_i)/a(num_i,num_j); while i>=1 s=0;

for k=i+1:num_i s=s+a(i,k)*x(k); end

x(i)=(b(i)-s)/a(i,i); i=i-1;

end case 1

x(1)=b(1)/a(1,1); j=2;

while j

s=s+a(j,k)*x(k); end

x(j)=(b(j)-s)/a(j,j); j=j+1; end otherwise

disp('error！请输入正确的参数！') end

运行结果：

⎛0. 0390

0. 0159 0. 0049

0. 0008A -1= 0. 0005

0. 0001 0 0. 0001 0. 0003⎝

0. 01600. 00580. 00360. 00730. 01760. 01290. 00140. 0012⎫

⎪

0. 03780. 01310. 00650. 01210. 00970. 00710. 00240. 0022⎪0. 01160. 03810. 00770. 00690. 00390. 00280. 00150. 0025⎪

⎪

0. 00200. 00640. 01810. 01050. 00330. 00240. 00240. 0058⎪0. 00110. 00360. 01010. 02430. 00700. 00510. 00480. 0035⎪⎪0. 00030. 00100. 00280. 00680. 04190. 03060. 00130. 0010⎪

⎪

0. 00020. 00070. 00210. 00500. 03060. 04680. 00100. 0007⎪0. 00020. 00080. 00230. 00480. 00140. 00100. 03820. 0033⎪

⎪

0. 00060. 00200. 00580. 00360. 00110. 00080. 00340. 0365⎪⎭

（4）求A 的行列式的计算程序及结果：

for i=1:num_i s=s*U(i,i); end

disp(['矩阵的行列式值为：',num2str(s)])

运行程序后结果与调用matlab 中det （）函数结果如下：

（1）求cholesky 分解程序及结果：

function l=chole(a) [num_i,num_j]=size(a); l=zeros(num_i); l(1,1)=sqrt(a(1,1)); for k=2:num_i

l(k,1)=a(k,1)/l(1,1); end j=2;

while j

l(j,j)=sqrt(a(j,j)-s1); i=j+1;

while i

s2=s2+l(i,k)*l(j,k); end

l(i,j)=(a(i,j)-s2)/l(j,j); i=i+1; end j=j+1; end

运行程序后的结果为：

000⎫⎛1. 4142

⎪

00⎪ 0. 70712. 1213

L = ⎪-0. 70711. 17852. 84800

⎪ 0. 70713. 0641-0. 74130. 7494⎪⎝⎭

（2）求解方程组程序及结果：

a=[2,1,-1,1;1,5,2,7;-1,2,10,1;1,7,1,11]; b=[13;-9;6;0]; L=chole(a); y=hd(1,L,b); x=hd(0,L',y); disp(x)

运行后的结果为:

=(26. 4146-65. 390212. 512238. 0732)

程序和结果如下：

A=[4,2,0,0;3,-2,1,0;0,2,5,3;0,0,-1,6]; B=[6;2;10;5]; [L,U]=sjfj(A);% y=hd(1,L,B); x=hd(0,U,y); disp('方程组的解为：')

disp(x)

function [l,u]=sjfj(a) num_i=size(a,1); i=2;

while i

a(i,i-1)=a(i,i-1)/a(i-1,i-1); a(i,i)=a(i,i)-a(i,i-1)*a(i-1,i); i=i+1; end l=tril(a); u=triu(a); for j=1:num_i; l(j,j)=1; end

运行程序后结果为：x =(1111)

编程求解矩阵A 的QR 分解：

（1）QR 分解计算程序：

function [Q,R]=hqr(A) [n,n]=size(A); Q=eye(n); for i=1:(n-1) E=eye(n-i+1); e1=E(:,1);

a=zeros(1,n-i+1)'; for j=1:(n-i+1) a(j)=A(j+i-1,i); end

av=sqrt(a'*a); w=a-av*e1;

h=E-(2/(w'*w))*(w*w'); q=eye(n); for k=i:n for j=i:n

q(k,j)=h(k-i+1,j-i+1); end end A=q*A; Q=q*Q; end R=A; Q=Q';

（2）输入矩阵A 后的计算结果：

⎛0. 31620. 70710. 2582-0. 5774⎫ ⎪ -0. 31620. 7071-0. 25820. 5774⎪

Q = ⎪-0. 632500. 77460

⎪ -0. 63250-0. 5164-0. 5774⎪⎝⎭⎛3. 1623-3. 1623-0. 4743-2. 0555⎫ ⎪

02. 8284-0. 35360. 3536 ⎪

R = ⎪001. 5492-1. 0328

⎪ 000-1. 1547⎪⎝⎭

（1）Jacobi 迭代法求解程序与结果：

a=0; b=0; c=0; while 1

a0=0.25*(7+b-c); b0=(1/8)*(-21-4*a-c); c0=(1/5)*(15+2*a-b);

A=[abs(a-a0),abs(b-b0),abs(c-c0)]; f=max(A); if f

disp('方程组的解为x') disp(x)

运行后的结果为：x =(0. 060-3. 1113. 647)

（2）Gauss-Seidel 迭代法求解的程序与结果：

a=0; b=0; c=0; while 1 a0=a; b0=b; c0=c;

a=(1/4)*(7+b-c); b=(1/8)*(-21-4*a-c); c=(1/5)*(15+2*a-b);

A=[abs(a-a0),abs(b-b0),abs(c-c0)]; f=max(A); if f

end end

disp('方程组的解为：') disp(x)

运行后的结果为：x =

(0. 061-3. 1113. 646)

（1）计算2x 3-5x -1=0在[1, 2]上的根的程序：

syms x f=2*x^3-5*x-1; df=diff(f,x); g=x-(f/df); x0=1; while 1

x1=subs(g,x0); dx=abs(x1-x0); if dx

disp(['Newton法的解为：x=',num2str(x1)]) break else x0=x1; end end x0=1; x1=1.1; while 1

x2=x1-(subs(f,x1)/(subs(f,x1)-subs(f,x0)))*(x1-x0); dx=abs(x2-x1); if dx

disp(['割线法的解为：x=',num2str(x2)]) break else x0=x1; x1=x2; end end

运行后结果为：Newton 法x =1. 673，割线法x =1. 673。

（2）计算e sin x =0在[-4, -3]上的根的程序：

syms x

f=exp(x)*sin(x); df=diff(f,x);

g=x-(f/df); x0=-2; while 1

x1=subs(g,x0); dx=abs(x1-x0); if dx

disp(['Newton法的解为：x=',num2str(x1)]) break else x0=x1; end end x0=-2; x1=-2.1; while 1

x2=x1-(subs(f,x1)/(subs(f,x1)-subs(f,x0)))*(x1-x0); dx=abs(x2-x1); if dx

disp(['割线法的解为：x=',num2str(x2)]) break else x0=x1; x1=x2; end end

运行后结果为：Newton 法x =-3. 1416，割线法x =-3. 1416。

syms x f=x*cos(x)-2; x0=-4; x2=-2; while 1 x1=0.5*(x0+x2); f0=subs(f,x0); f1=subs(f,x1); c=f0*f1; if c

l=abs(x2-x0); if l

disp(['二分法的解为x=',num2str(x0)]) break end end

运行后的结果为x =-2. 499。

（1）Lagrange 插值程序：

function yv=lagran(xd,yd,xv) syms x

num=length(xd);%计算节点的个数 f=0;

for k=1:num%创建lagrange 插值函数 index=setdiff(1:num,k);

f=f+yd(k)*prod((x-xd(index))./(xd(k)-xd(index))); end

yv=subs(f,xv);%计算xv 点插值函数值

（2）绘制函数图程序： syms x

f=1./(1+x.^2); xd2=-5:2:5;

yd2=subs(f,xd2); x0=-5:0.1:5; y=subs(f,x0); plot(x0,y) hold on

y2=lagran(xd2,yd2,x0); plot(x0,y2)

（3）运行程序后得到的步长分别为2，1，

的插值函数图与原图形的比较图如下： 2

（1）三次样条插值程序：

function yv=csi(xd,yd,xv,h,mf,ml) num_xd=length(xd); n=num_xd-1; a_m=1:num_xd-2; for i=1:num_xd-2

a_m(i)=2; end

a_np=1:num_xd-3; for i=1:num_xd-3 a_np(i)=0.5; end

A=diag(a_m)+diag(a_np,-1)+diag(a_np,1); g=1:n-1; for k=1:n-1

g(k)=3*((0.5/h)*(yd(k+2)-yd(k+1))+(0.5/h)*(yd(k+1)-yd(k))); end b=1:n-1;

b(1)=g(1)-0.5*mf; b(n-1)=g(n-1)-0.5*ml; for k=2:n-2 b(k)=g(k); end b=b.'; M=A\b; m=1:n+1; m(1)=mf; m(n+1)=ml; for i=2:n m(i)=M(i-1); end

num_xv=length(xv); yv=1:num_xv; for i=1:num_xv for j=1:num_xd if xv(i)

yv1=(h+2*(xv(i)-xd(k1)))*(xv(i)-xd(k2))^2*(yd(k1)/h^3); yv2=(h-2*(xv(i)-xd(k2)))*(xv(i)-xd(k1))^2*(yd(k2)/h^3); yv3=(xv(i)-xd(k1))*(xv(i)-xd(k2))^2*(m(k1)/h^2); yv4=(xv(i)-xd(k2))*(xv(i)-xd(k1))^2*(m(k2)/h^2); yv(i)=yv1+yv2+yv3+yv4; end

（2）绘图程序（改变程序中的h 值即可改变步长）：

clear all clc

h=0.5;

xd=-5:h:5;

yd=1./(1+xd.^2);

xv=-5:0.1:5;

y=1./(1+xv.^2);

mf=0.014793;

ml=-0.014793;

yv=csi(xd,yd,xv,h,mf,ml);

plot(xv,y)

hold on

plot(xv,yv)

（3）步长为2，1，1时原函数图与插值函数图的比较：

（1）复化梯形公式计算积分程序：

function t=trape(f,a,b,n)

h=(b-a)/n;

index=(a+h):h:(b-h);

t1=sum(subs(f,index));

t=((b-a)/(2*n))*(subs(f,a)+2*t1+subs(f,b));

（2）复化Simpson 公式计算程序：

function s=simpson(f,a,b,n)

h=(b-a)/n;

index1=(a+0.5*h):h:(b-0.5*h);

index2=(a+h):h:(b-h);

s1=sum(subs(f,index1));

s2=sum(subs(f,index2));

s=((b-a)/(6*n))*(subs(f,a)+4*s1+2*s2+subs(f,b));

2（3）区间数不同时的计算结果：（其中⎰-2x 2sinxdx =0）

（1）Gauss 积分法计算积分程序：

function s=gauss(f,a,b,n)

switch n%Gauss

case 2

x=[-1/sqrt(3),1/sqrt(3)];

w=[1,1];

case 3

x=[-0.7745966692,0,0.7745966692];

w=[0.5555555556,0.8888888889,0.5555555556];

case 4

x=[-0.8611363116,-0.3399810436,0.3399810436,0.8611363116]; w=[0.3478548451,0.6521451549,0.6521451549,0.3478548451]; case 5

x=[-0.9061798459,-0.5384693101,0,0.5384693101,0.9061798459];

w=[0.2369268851,0.4786286705,0.5688888889,0.4786286705,0.2369268851]; end

T=(a+b)/2+(b-a)/2*x;

s=(b-a)/2*sum(w.*subs(f,T));

（2）当点数为2，3，5时的计算结果：

当点数为2时计算结果：s =0. 4746

当点数为3时计算结果：s =0. 4672

当点数为5时计算结果：s =0. 4674

由计算结果可知点数越多，计算越精确。

（1）Euler 法计算程序：

function [T,U]=euler(f,a,b,ul,h)

syms t u

n=(b-a)/h;

U=1:(n+1);

T=a:h:b;

U(1)=ul;

for i=1:

U(i+1)=U(i)+h*subs(f,{t,u},{T(i),U(i)});

end

（2）改进的Euler 法计算程序：

function [T,U]=imeuler(f,a,b,ul,h)

syms t u

n=(b-a)/h;

T=a:h:b;

U=1:(n+1);

U(1)=ul;

for i=1:n

m=subs(f,{t,u},{T(i),U(i)});

U(i+1)=U(i)+0.5*h*(m+subs(f,{t,u},{T(i+1),U(i)+h*m}));

end

（3）Runge-Kutta 法计算程序：

function [T,U]=ruku(f,a,b,ul,h)

syms t u

n=(b-a)/h;

T=a:h:b;

U=1:(n+1);

U(1)=ul;

for i=1:n

k1=subs(f,{t,u},{T(i),U(i)});

k2=subs(f,{t,u},{T(i)+0.5*h,U(i)+0.5*h*k1}); k3=subs(f,{t,u},{T(i)+0.5*h,U(i)+0.5*h*k2}); k4=subs(f,{t,u},{T(i)+h,U(i)+h*k3});

U(i+1)=U(i)+(h/6)*(k1+2*k2+2*k3+k4);

end

（4）当取步长分别为0.1，0.05，0.01时的各种方法的结果比较图：

（1）从大到小的顺序的计算程序：

function y=snd(n) format long y=0; if n

disp('请输入大于1的数！') else s=0; i=2;

while i

s=single(s+(1/(i^2-1))); i=i+1; end y=s; end

（2）从小到大的顺序的计算程序：

function y=snx(n) format long y=0; if n

disp('请输入大于1的数！') else s=0; i=n; while 1

s=single(s+(1/(i^2-1))); i=i-1; if i==1 break end end

y=s; end

（3）按两种顺序分别计算① S 102的计算结果：

② S 104的计算结果：

③ S 106的计算结果：

计算时的有效位数为七位数。

并指出有效位数（编制程序时用单精度）

① 秦九昭算法计算程序：

function y=qjz(a,x) j=3;

i=size(a,2); switch i case 1 y=a(1); case 2

y=a(1)*x+a(2); otherwise

p=a(1)*x+a(2); while j

② 计算

计算结果如下：

在点23的值。

当x =23时f (x )=86652。

① Gauss 法计算程序和结果：程序：

s=s-x(k)*a(j,k); k=k+1; end

x(j)=s/a(j,j); j=j-1; end disp(x)

function [x,y]=gauss(a,b) num_i=size(a,1); j=1;

while j

while i

b(i,:)=b(i,:)-r*b(j,:); i=i+1; end j=j+1; end x=a; y=b;

运行的结果为：

x =(-0. 2890. 345-0. 713-0. 221-0. 4300. 154-0. 0580. 2010. 290)。

② 列主元消去法计算程序和结果：

A(1,:)=[31,-13,0,0,0,-10,0,0,0]; A(2,:)=[-13,35,-9,0,-11,0,0,0,0]; A(3,:)=[0,-9,31,-10,0,0,0,0,0]; A(4,:)=[0,0,-10,79,-30,0,0,0,-9]; A(5,:)=[0,0,0,-30,57,-7,0,-5,0]; A(6,:)=[0,0,0,0,-7,47,-30,0,0]; A(7,:)=[0,0,0,0,0,-30,41,0,0]; A(8,:)=[0,0,0,0,-5,0,0,27,-2]; A(9,:)=[0,0,0,-9,0,0,0,-2,29]; B=[-15;27;-23;0;-20;12;-7;7;10]; [a,b]=lzy(A,B); j=size(a,2); while j>=1 k=j+1; s=b(j); while k

s=s-x(k)*a(j,k); k=k+1; end

x(j)=s/a(j,j); j=j-1; end disp(x)

function [a1,b1]=lzy(a,b) [num_i,num_j]=size(a); ab=zeros(num_i,num_j+1); for k=1:num_j ab(:,k)=a(:,k); end

ab(:,num_j+1)=b(:,1); j=1;

while j

[max,max_i]=searmax(j,j,ab); i_nu=ab(j,:); ab(j,:)=ab(max_i,:); ab(max_i,:)=i_nu; m=j+1;

while m

ab(m,n)=ab(m,n)-(ab(m,j)/max)*ab(j,n); end m=m+1; end j=j+1; end

a1=zeros(num_i,num_j); for k=1:num_i a1(:,k)=ab(:,k); end

b1=ab(:,num_i+1);

function [b,c]=searmax(i,j,a) num_i=size(a,1); k=i;

m=abs(a(k,j)); c=k;

while k

if m>=abs(a(k,j)) continue else

m=abs(a(k,j)); c=k; end end b=a(c,j);

运行的结果为：

x =

(-0. 2370. 477-0. 767-0. 078-0. 3080. 146-0. 1710. 2850. 345)

（1）LU 分解的计算程序及结果：

function [l,u]=lufz(a) [num_i,num_j]=size(a);

l=eye(num_i); u=eye(num_i); for j=1:num_j u(1,j)=a(1,j); l(j,1)=a(j,1)/u(1,1); end i=2;

while i

while j

s=s+l(i,k)*u(k,j); end

u(i,j)=a(i,j)-s; s=0; for k=1:i-1

s=s+l(j+1,k)*u(k,i); end

l(j+1,i)=(a(j+1,i)-s)/u(i,i); j=j+1; end s=0; for k=1:i-1

s=s+l(i,k)*u(k,num_i); end

u(i,num_i)=a(i,num_i)-s; i=i+1; end

输入要求解的矩阵后运行的结果为：

⎛ 1 -0. 4191 0-0. 3051 00-0. 3541L = 000-0. 3981 0000-0. 1571 00000-0. 6551 0000-0. 112-0. 018-0. 0251 ⎝000-0. 119-0. 083-0. 014-0. 020-0. 092⎫⎪⎪⎪⎪

⎪

⎪⎪

⎪⎪⎪⎪⎭

29. 548-90-11 28. 259-10-3. 350

75. 461-31. 186

U = 44. 602

⎝

⎪⎪

-1. 277000⎪

⎪

-0. 45200-9⎪-7. 1800-5-3. 578⎪⎪45. 8732-30-0. 785-0. 562⎪

⎪

21. 381-0. 513-0. 367⎪

26. 413-2. 420⎪

⎪27. 390⎪⎭

-4. 194000

（2）带列主元的LU 分解计算程序及结果

由于Matlab 中自带带列主元的LU 分解函数，故计算程序如下：

运行结果如下：

⎛1⎫ ⎪

1 -0. 419⎪

0⎪-0. 3051 ⎪

0-0. 3541 0⎪

⎪

L = 000-0. 3981 ⎪

0⎪000-0. 1571 ⎪

00000-0. 6551 ⎪ 0⎪000-0. 112-0. 018-0. 0251 ⎪ 0⎪00-0. 119-0. 083-0. 014-0. 020-0. 0921⎝⎭

29. 548-90-11 28. 259-10-3. 350

75. 461-31. 186

U = 44. 602

⎪⎪

-1. 277000⎪

⎪

-0. 45200-9⎪-7. 1800-5-3. 578⎪⎪45. 8732-30-0. 785-0. 562⎪

⎪

21. 381-0. 513-0. 367⎪

-4. 194000

⎝⎛ 100000000⎫

010000000⎪

⎪ 001000000⎪

000100000⎪

⎪P = 0

00010000⎪ 000001000⎪为单位阵。 ⎪

000000100⎪

⎪ 000000010⎪

⎝

000000001⎪

⎪⎭

（3）求A 的逆矩阵的计算程序及结果：

A(1,:)=[31,-13,0,0,0,-10,0,0,0]; A(2,:)=[-13,35,-9,0,-11,0,0,0,0]; A(3,:)=[0,-9,31,-10,0,0,0,0,0]; A(4,:)=[0,0,-10,79,-30,0,0,0,-9]; A(5,:)=[0,0,0,-30,57,-7,0,-5,0]; A(6,:)=[0,0,0,0,-7,47,-30,0,0]; A(7,:)=[0,0,0,0,0,-30,41,0,0]; A(8,:)=[0,0,0,0,-5,0,0,27,-2]; A(9,:)=[0,0,0,-9,0,0,0,-2,29]; [num_i,num_j]=size(A); A_n=eye(num_i); E=eye(num_i); [L,U]=lufz(A); for i=1:num_i y=hd(1,L,E(:,i)); A_n(:,i)=hd(0,U,y); end disp(A_n)

function [l,u]=lufz(a) [num_i,num_j]=size(a);

26. 413-2. 420⎪

27. 390⎪⎪⎭

l=eye(num_i); u=eye(num_i); for j=1:num_j u(1,j)=a(1,j); l(j,1)=a(j,1)/u(1,1); end i=2;

while i

while j

s=s+l(i,k)*u(k,j); end

u(i,j)=a(i,j)-s; s=0; for k=1:i-1

s=s+l(j+1,k)*u(k,i); end

l(j+1,i)=(a(j+1,i)-s)/u(i,i); j=j+1; end s=0; for k=1:i-1

s=s+l(i,k)*u(k,num_i); end

u(i,num_i)=a(i,num_i)-s; i=i+1; end

function x=hd(f,a,b) [num_i,num_j]=size(a); x=zeros(num_i,1); switch f case 0 i=num_i-1;

x(num_i)=b(num_i)/a(num_i,num_j); while i>=1 s=0;

for k=i+1:num_i s=s+a(i,k)*x(k); end

x(i)=(b(i)-s)/a(i,i); i=i-1;

end case 1

x(1)=b(1)/a(1,1); j=2;

while j

s=s+a(j,k)*x(k); end

x(j)=(b(j)-s)/a(j,j); j=j+1; end otherwise

disp('error！请输入正确的参数！') end

运行结果：

⎛0. 0390

0. 0159 0. 0049

0. 0008A -1= 0. 0005

0. 0001 0 0. 0001 0. 0003⎝

0. 01600. 00580. 00360. 00730. 01760. 01290. 00140. 0012⎫

⎪

0. 03780. 01310. 00650. 01210. 00970. 00710. 00240. 0022⎪0. 01160. 03810. 00770. 00690. 00390. 00280. 00150. 0025⎪

⎪

0. 00200. 00640. 01810. 01050. 00330. 00240. 00240. 0058⎪0. 00110. 00360. 01010. 02430. 00700. 00510. 00480. 0035⎪⎪0. 00030. 00100. 00280. 00680. 04190. 03060. 00130. 0010⎪

⎪

0. 00020. 00070. 00210. 00500. 03060. 04680. 00100. 0007⎪0. 00020. 00080. 00230. 00480. 00140. 00100. 03820. 0033⎪

⎪

0. 00060. 00200. 00580. 00360. 00110. 00080. 00340. 0365⎪⎭

（4）求A 的行列式的计算程序及结果：

for i=1:num_i s=s*U(i,i); end

disp(['矩阵的行列式值为：',num2str(s)])

运行程序后结果与调用matlab 中det （）函数结果如下：

（1）求cholesky 分解程序及结果：

function l=chole(a) [num_i,num_j]=size(a); l=zeros(num_i); l(1,1)=sqrt(a(1,1)); for k=2:num_i

l(k,1)=a(k,1)/l(1,1); end j=2;

while j

l(j,j)=sqrt(a(j,j)-s1); i=j+1;

while i

s2=s2+l(i,k)*l(j,k); end

l(i,j)=(a(i,j)-s2)/l(j,j); i=i+1; end j=j+1; end

运行程序后的结果为：

000⎫⎛1. 4142

⎪

00⎪ 0. 70712. 1213

L = ⎪-0. 70711. 17852. 84800

⎪ 0. 70713. 0641-0. 74130. 7494⎪⎝⎭

（2）求解方程组程序及结果：

a=[2,1,-1,1;1,5,2,7;-1,2,10,1;1,7,1,11]; b=[13;-9;6;0]; L=chole(a); y=hd(1,L,b); x=hd(0,L',y); disp(x)

运行后的结果为:

=(26. 4146-65. 390212. 512238. 0732)

程序和结果如下：

A=[4,2,0,0;3,-2,1,0;0,2,5,3;0,0,-1,6]; B=[6;2;10;5]; [L,U]=sjfj(A);% y=hd(1,L,B); x=hd(0,U,y); disp('方程组的解为：')

disp(x)

function [l,u]=sjfj(a) num_i=size(a,1); i=2;

while i

a(i,i-1)=a(i,i-1)/a(i-1,i-1); a(i,i)=a(i,i)-a(i,i-1)*a(i-1,i); i=i+1; end l=tril(a); u=triu(a); for j=1:num_i; l(j,j)=1; end

运行程序后结果为：x =(1111)

编程求解矩阵A 的QR 分解：

（1）QR 分解计算程序：

function [Q,R]=hqr(A) [n,n]=size(A); Q=eye(n); for i=1:(n-1) E=eye(n-i+1); e1=E(:,1);

a=zeros(1,n-i+1)'; for j=1:(n-i+1) a(j)=A(j+i-1,i); end

av=sqrt(a'*a); w=a-av*e1;

h=E-(2/(w'*w))*(w*w'); q=eye(n); for k=i:n for j=i:n

q(k,j)=h(k-i+1,j-i+1); end end A=q*A; Q=q*Q; end R=A; Q=Q';

（2）输入矩阵A 后的计算结果：

⎛0. 31620. 70710. 2582-0. 5774⎫ ⎪ -0. 31620. 7071-0. 25820. 5774⎪

Q = ⎪-0. 632500. 77460

⎪ -0. 63250-0. 5164-0. 5774⎪⎝⎭⎛3. 1623-3. 1623-0. 4743-2. 0555⎫ ⎪

02. 8284-0. 35360. 3536 ⎪

R = ⎪001. 5492-1. 0328

⎪ 000-1. 1547⎪⎝⎭

（1）Jacobi 迭代法求解程序与结果：

a=0; b=0; c=0; while 1

a0=0.25*(7+b-c); b0=(1/8)*(-21-4*a-c); c0=(1/5)*(15+2*a-b);

A=[abs(a-a0),abs(b-b0),abs(c-c0)]; f=max(A); if f

disp('方程组的解为x') disp(x)

运行后的结果为：x =(0. 060-3. 1113. 647)

（2）Gauss-Seidel 迭代法求解的程序与结果：

a=0; b=0; c=0; while 1 a0=a; b0=b; c0=c;

a=(1/4)*(7+b-c); b=(1/8)*(-21-4*a-c); c=(1/5)*(15+2*a-b);

A=[abs(a-a0),abs(b-b0),abs(c-c0)]; f=max(A); if f

end end

disp('方程组的解为：') disp(x)

运行后的结果为：x =

(0. 061-3. 1113. 646)

（1）计算2x 3-5x -1=0在[1, 2]上的根的程序：

syms x f=2*x^3-5*x-1; df=diff(f,x); g=x-(f/df); x0=1; while 1

x1=subs(g,x0); dx=abs(x1-x0); if dx

disp(['Newton法的解为：x=',num2str(x1)]) break else x0=x1; end end x0=1; x1=1.1; while 1

x2=x1-(subs(f,x1)/(subs(f,x1)-subs(f,x0)))*(x1-x0); dx=abs(x2-x1); if dx

disp(['割线法的解为：x=',num2str(x2)]) break else x0=x1; x1=x2; end end

运行后结果为：Newton 法x =1. 673，割线法x =1. 673。

（2）计算e sin x =0在[-4, -3]上的根的程序：

syms x

f=exp(x)*sin(x); df=diff(f,x);

g=x-(f/df); x0=-2; while 1

x1=subs(g,x0); dx=abs(x1-x0); if dx

disp(['Newton法的解为：x=',num2str(x1)]) break else x0=x1; end end x0=-2; x1=-2.1; while 1

x2=x1-(subs(f,x1)/(subs(f,x1)-subs(f,x0)))*(x1-x0); dx=abs(x2-x1); if dx

disp(['割线法的解为：x=',num2str(x2)]) break else x0=x1; x1=x2; end end

运行后结果为：Newton 法x =-3. 1416，割线法x =-3. 1416。

syms x f=x*cos(x)-2; x0=-4; x2=-2; while 1 x1=0.5*(x0+x2); f0=subs(f,x0); f1=subs(f,x1); c=f0*f1; if c

l=abs(x2-x0); if l

disp(['二分法的解为x=',num2str(x0)]) break end end

运行后的结果为x =-2. 499。

（1）Lagrange 插值程序：

function yv=lagran(xd,yd,xv) syms x

num=length(xd);%计算节点的个数 f=0;

for k=1:num%创建lagrange 插值函数 index=setdiff(1:num,k);

f=f+yd(k)*prod((x-xd(index))./(xd(k)-xd(index))); end

yv=subs(f,xv);%计算xv 点插值函数值

（2）绘制函数图程序： syms x

f=1./(1+x.^2); xd2=-5:2:5;

yd2=subs(f,xd2); x0=-5:0.1:5; y=subs(f,x0); plot(x0,y) hold on

y2=lagran(xd2,yd2,x0); plot(x0,y2)

（3）运行程序后得到的步长分别为2，1，

的插值函数图与原图形的比较图如下： 2

（1）三次样条插值程序：

function yv=csi(xd,yd,xv,h,mf,ml) num_xd=length(xd); n=num_xd-1; a_m=1:num_xd-2; for i=1:num_xd-2

a_m(i)=2; end

a_np=1:num_xd-3; for i=1:num_xd-3 a_np(i)=0.5; end

A=diag(a_m)+diag(a_np,-1)+diag(a_np,1); g=1:n-1; for k=1:n-1

g(k)=3*((0.5/h)*(yd(k+2)-yd(k+1))+(0.5/h)*(yd(k+1)-yd(k))); end b=1:n-1;

b(1)=g(1)-0.5*mf; b(n-1)=g(n-1)-0.5*ml; for k=2:n-2 b(k)=g(k); end b=b.'; M=A\b; m=1:n+1; m(1)=mf; m(n+1)=ml; for i=2:n m(i)=M(i-1); end

num_xv=length(xv); yv=1:num_xv; for i=1:num_xv for j=1:num_xd if xv(i)

（2）绘图程序（改变程序中的h 值即可改变步长）：

clear all clc

h=0.5;

xd=-5:h:5;

yd=1./(1+xd.^2);

xv=-5:0.1:5;

y=1./(1+xv.^2);

mf=0.014793;

ml=-0.014793;

yv=csi(xd,yd,xv,h,mf,ml);

plot(xv,y)

hold on

plot(xv,yv)

（3）步长为2，1，1时原函数图与插值函数图的比较：

（1）复化梯形公式计算积分程序：

function t=trape(f,a,b,n)

h=(b-a)/n;

index=(a+h):h:(b-h);

t1=sum(subs(f,index));

t=((b-a)/(2*n))*(subs(f,a)+2*t1+subs(f,b));

（2）复化Simpson 公式计算程序：

function s=simpson(f,a,b,n)

h=(b-a)/n;

index1=(a+0.5*h):h:(b-0.5*h);

index2=(a+h):h:(b-h);

s1=sum(subs(f,index1));

s2=sum(subs(f,index2));

s=((b-a)/(6*n))*(subs(f,a)+4*s1+2*s2+subs(f,b));

2（3）区间数不同时的计算结果：（其中⎰-2x 2sinxdx =0）

（1）Gauss 积分法计算积分程序：

function s=gauss(f,a,b,n)

switch n%Gauss

case 2

x=[-1/sqrt(3),1/sqrt(3)];

w=[1,1];

case 3

x=[-0.7745966692,0,0.7745966692];

w=[0.5555555556,0.8888888889,0.5555555556];

case 4

x=[-0.8611363116,-0.3399810436,0.3399810436,0.8611363116]; w=[0.3478548451,0.6521451549,0.6521451549,0.3478548451]; case 5

x=[-0.9061798459,-0.5384693101,0,0.5384693101,0.9061798459];

w=[0.2369268851,0.4786286705,0.5688888889,0.4786286705,0.2369268851]; end

T=(a+b)/2+(b-a)/2*x;

s=(b-a)/2*sum(w.*subs(f,T));

（2）当点数为2，3，5时的计算结果：

当点数为2时计算结果：s =0. 4746

当点数为3时计算结果：s =0. 4672

当点数为5时计算结果：s =0. 4674

由计算结果可知点数越多，计算越精确。

（1）Euler 法计算程序：

function [T,U]=euler(f,a,b,ul,h)

syms t u

n=(b-a)/h;

U=1:(n+1);

T=a:h:b;

U(1)=ul;

for i=1:

U(i+1)=U(i)+h*subs(f,{t,u},{T(i),U(i)});

end

（2）改进的Euler 法计算程序：

function [T,U]=imeuler(f,a,b,ul,h)

syms t u

n=(b-a)/h;

T=a:h:b;

U=1:(n+1);

U(1)=ul;

for i=1:n

m=subs(f,{t,u},{T(i),U(i)});

U(i+1)=U(i)+0.5*h*(m+subs(f,{t,u},{T(i+1),U(i)+h*m}));

end

（3）Runge-Kutta 法计算程序：

function [T,U]=ruku(f,a,b,ul,h)

syms t u

n=(b-a)/h;

T=a:h:b;

U=1:(n+1);

U(1)=ul;

for i=1:n

k1=subs(f,{t,u},{T(i),U(i)});

k2=subs(f,{t,u},{T(i)+0.5*h,U(i)+0.5*h*k1}); k3=subs(f,{t,u},{T(i)+0.5*h,U(i)+0.5*h*k2}); k4=subs(f,{t,u},{T(i)+h,U(i)+h*k3});

U(i+1)=U(i)+(h/6)*(k1+2*k2+2*k3+k4);

end

（4）当取步长分别为0.1，0.05，0.01时的各种方法的结果比较图：

大连理工大学-矩阵大作业

相关内容

热门内容

标签