北京科技大学应用计算方法作业与答案

一、第一次作业

（一）2-6计算下列向量的1-范数、∞-范数、2-范数。

（1）x=(12,-4,-6,2)T

>> A=[12,-4,-6,2] A =

12 -4 -6 2 >> norm(A,1) ans = ^

24

>> norm(A,inf) ans = 12

>> norm(A,2) ans =

(2) x=(1,3,-4)T

、

>> A=[1,3,-4] A =

1 3 -4 >> norm(A,1) ans = 8

>> norm(A,inf) ans = …

4

>> norm(A,2) ans =

（二）2-9 计算下列矩阵的行范数、列范数、谱范数、F 范数。

（1）???

?

? ??--=11211

1113A >> A=[3,-1,1;1,1,1;2,1,-1]

A = ~

3 -1 1 1 1 1

2 1 -1 >> norm(A,1) ans = 6

>> norm(A,inf) ans = …

5

>> norm(A,2) ans =

>> norm(A,'fro') ans = （2）R a a a A ∈?

??

?

??-=,00 .

>> A=[0,1;-1,0] A =

0 1 -1 0 >> norm(A,1) ans = 1

>> norm(A,inf) [

ans = 1

>> norm(A,2) ans = 1

>> norm(A,'fro') ans =

>

二、第二次作业

用牛顿迭代法求方程0133=--x x 在20=x 附近的根。要求：给成程序和运行结

果.

1、牛顿法的基本原理

在求解非线性方程0)(=x f 时，它的困难在于)(x f 是非线性函数，为克服这一困难，考虑它的线性展开。设当前点为k x ，在k x 处的Taylor 展开式为

))(()()('

k k k x x x f x f x f -+≈ )1.2(

令0)(=x f ，解其方程得到

),1,0(,)

()

('1 =-

=+k x f x f x x k k k k )2.2(

)2.2(式为牛顿迭代公式，用牛顿迭代公式求方程0)(=x f 根的方法称为牛顿迭代法。此即

牛顿迭代法的设计原理。

2、matlab 程序代码

function root=NewtonRoot(f,a,b,eps) if (nargin==3) eps=; end

f1=subs(sym(f),findsym(sym(f)),a); f2=subs(sym(f),findsym(sym(f)),b); if (f1==0) root=a; end

if (f2==0) roor=b; end

if (f1*f2>0)

disp('两端点函数值乘积大于0！') return ; else

tol=1;

fun=diff(sym(f));

fa=subs(sym(f),findsym(sym(f)),a); fb=subs(sym(f),findsym(sym(f)),b);

dfa=subs(sym(fun),findsym(sym(fun)),a); dfb=subs(sym(fun),findsym(sym(fun)),b); if (dfa>dfb)

root=a-fa/dfa; else

root=b-fb/dfb; end

while(tol>eps)

r1=root;

fx=subs(sym(f),findsym(sym(f)),r1);

dfx=subs(sym(fun),findsym(sym(fun)),r1);

root=r1-fx/dfx;

tol=abs(root-r1);

end

end

3、运行结果截图

结论：通过计算可以看出0

1

3

3=

-

-x

x在2

=

x附近的根为.

三、第三次作业

编写高斯顺序消元法求解下面方程组的程序并计算结果。

?

?

?

?

?

=

+

-

-

=

-

+

-

=

-

-

2.3

5

2

5.8

2

10

2.6

10

3

2

1

3

2

1

3

2

1

x

x

x

x

x

x

x

x

x

)1.3(

1、高斯顺序消元法的设计原理

高斯顺序消元法的基本思想是将线性方程组

?

?

?

?

?

?

?

=

+

+

+

=

+

+

+

=

+

+

+

n

n

nn

n

n

n

n

n

n

b

x

a

x

a

x

a

b

x

a

x

a

x

a

b

x

a

x

a

x

a

2

2

1

1

2

2

2

22

1

21

1

1

2

12

1

11

)2.3(

通过消元，逐步转化为等价的上（或下）三角形方程组，然后用回代法求解。

2、matlab程序代码

function [x,XA]=GaussXQByOrder(A,b)

%高斯顺序消元法

N = size(A);

n = N(1);

for i=1:(n-1)

for j=(i+1):n

if(A(i,i)==0)

disp('对角元素为0！'); %防止对角元素为0

return;

end

l = A(j,i);

m = A(i,i);

A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m; %消元方程

b(j)=b(j)-l*b(i)/m;

end

end

x=SolveUpTriangle(A,b); %通用的求上三角系数矩阵线性方程组的函数XA = A; %消元后的系数矩阵

function x=SolveUpTriangle(A,b)

N=size(A);

n=N(1);

for i=n:-1:1

if(i

s=A(i,(i+1):n)*x((i+1):n,1);

else

s=0;

end

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

end

3、运行结果截图

结论：高斯顺序消元法求解出)1.3(的结果为2200

.1

,

1800

.1

,

8600

.0

3

2

1

=

=

=x

x

x。

四、第四次作业

编写Jacobi迭代法和Seidel迭代法求解方程组的程序，并计算出结果。

?

?

?

?

?

=

+

-

-

=

-

+

-

=

-

-

2.3

5

2

5.8

2

10

2.6

10

3

2

1

3

2

1

3

2

1

x

x

x

x

x

x

x

x

x

)1.4(精度要求：3

)

(

)1

(10-

∞

+<

-k

k x

x

（一）用Jacobi求解题设方程组

1、Jacobi迭代原理

设有一个n元线性方程组

??

????

?=+++=+++=+++n

n nn n n n n n n b x a x a x a b x a x a x a b x a x a x a 22112

222212********* )2.4( 它的矩阵形式为B AX =，如果n n ij a A ?=)(非奇异，且n i a ii ,2,1,0=≠。由上式方程组可以得到

),,2,1(),(11

1

n i x a b a x n

j j j ij i ii i =-=∑≠= )3.4(

而其相应的迭代公式

),,2,1(),(11

1

)(n i x a b a x n j j k j ij i ii i =-=∑≠= )4.4(

此迭代公式即为Jacobi 迭代。

2、Jacobi 迭代法matlab 程序代码

function [x,n]=jacobi(A,b,x0,eps,varargin)

% 求解线性方程组的迭代法其中 % A 为方程组的系数矩阵 % b 为方程组的右端项

% eps 为精度要求，默认值为1e-5 % varargin 为最大迭代次数，值100 % x 为方程组的解 % n 为迭代次数 if nargin==3 eps=; M=200;

elseif nargin<3 error return

elseif nargin==5 M=varargin{1}; end

D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=D\(L+U); f=D\b; x=B*x0+f; n=1;

while norm(x-x0)>=eps x0=x; x=B*x0+f;

n=n+1;

if(n>=M)

disp('Warning:迭代次数太多，可能不收敛！');

return;

end

end

3、运行结果截图

结论：)1.4(式方程组的解为2200

.1

,

1800

.1

,

8600

.0

3

2

1

=

=

=x

x

x，需要迭代12步。（二）用Seidel求解题设方程组

1、Seidel迭代原理

在Jacobi迭代计算过程中可看出，Jacobi方法计算)1

(+

k

i

x时，并未用到已算出的

)1

(

1

)1

(

2

)1

(

1

,

,+

-

+

+k

i

k

k x

x

x ，这时想到，如果Jacobi迭代收敛，)1

,

,1

()1

(-

=

+i

j

x k

j

比

)1

,

,1

()(-

=i

j

x k

j

更接近方程组的解，若能在迭代过程中尽快用新的信息)1

(+

k

j

x去替换)

(k

j

x，则可望收敛更快。由此，可将Jacobi迭代公式改写为

,2,1,0

,

,

,2,1

),

(

1

1

)

(

1

1

)1

(

)1

(=

=

-

-

=∑

∑

+

=

-

=

+

+k

n

i

x

a

x

a

b

a

x

n

i

j

k

j

ij

i

j

k

j

ij

i

ii

k

i

)5.4(

此式即为Seidel迭代法。

2、Seidel迭代法matlab程序代码

function[x,n]=gauseidel(A,b,x0,eps,M)

if nargin==3

eps=;

M=200;

elseif nargin==4 M=200;

elseif nargin<3 error return ; end

D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); G=(D-L)\U; f=(D-L)\b; x=G*x0+f; n=1;

while norm(x-x0)>=eps x0=x; x=G*x0+f; n=n+1; if (n>=M)

disp('Warning:迭代次数太多，可能不收敛！'); return ; end end

3、运行结果截图

结论：)1.4(式方程组的解为2200.1,1800.1,8600.0321===x x x ，需要迭代5步。

五、第五次作业

用归一化算法（归一化幂法6-1）求矩阵A 的最大模特征值和特征向量，其中

T V A )1,1,1(,322244461)

0(=????

? ??=

1、幂法的基本思想

任取一个非零初始向量n R v ∈0且00≠v ，由矩阵A 的乘幂构造一迭代序列为

??

????

?======-010

21201v

A Av v v A Av v Av v k k k

)1.5( 假设矩阵A 有n 个线性无关的特征向量),,2,1(n k x k =，于是给定的初始向量0v 可以用这组特征向量线性表示，即∑==+++=n

i i

i n n x

a x a x a x a v 1

22110 并设0≠i a 。把0v 代入

迭代序列的第一条式子，得

∑==+++=+++==n

i i i i n n n n n x a x a x a x a x a x a x a A Av v 1

222111221101)(λλλλ

同理可得

),,2,1(,1

22

211

1n k x a x a x a x a Av v n

i i k i i n k n

n k k

k k ==+++==∑=λλλλ

2、matlab 程序代码

function [l,v,s]=pmethod(A,x0,eps) if nargin==2 eps = ; end

v = x0; %v 为主特征向量 M = 5000; %迭代步数限制 m = 0; l = 0; for (k=1:M) y = A*v;

m = max(y); %m 为按模最大的分量 v = y/m;

if (abs(m - l)

l = m; %到所需精度，退出，l 为主特征值 s = k; %s 为迭代步数 return ; else

if (k==M)

disp('迭代步数太多，收敛速度太慢！'); l = m; s = M; else

l = m; end end end

3、运行结果截图

结论：由输出结果可知，经过15步迭代，求得矩阵特征值为，对应的迭代向量为（,，）T

六、第六次作业

已知)(x f 的函数值和导数值如下：2)2(,1)1(,3)3(,2)2(,1)1(''=====f f f f f 求次数小于等于4的多项式)(x P ，使2)2(,1)1(,3)3(,2)2(,1)1(''=====f f f f f 并给出余项公式。

解：首先构造拉格朗日插值函数

3)

23()13()

2)(1(2)32()12()3)(1(1)31()21()3)(2(2?-?---+?-?---+?-?---=

x x x x x x P

化简得x P =2

构造函数)3)(2)(1)((4---++=x x x b ax x P

对上式求导得1611221218342

2

3

'

4+-++--+=a b ax bx ax bx ax P 根据导数信息可以解出1,1=-=b a 所以)3)(2)(1)(1(4---+-+=x x x x x P

余项)3()2()1)(()()()(2

2

44---=-=x x x x k x P x f x R 构造辅助函数)3()2()1)(()()()(2

2

4-----=t t t x k t P t f t ? 其中)

3()2()1()

()(224---=

x x x x R x k

然后反复利用Rolle 定理最好可以得到

)3()2()1(!

5)

()(2254---=x x x f x R ξ

七、第七次作业

已知某地区在不同月份的平均日照时间的观测数据如下表，试分析日照时间的变化规律。

表1：日照时间

解：首先在matlab 中输入语句 >> x=[1:12];

>> y=[ ]; >> xx=::;

>> yy=spline(x,y,xx); >> plot(x,y,'o',xx,yy)

图1

通过图像可以看出日照随时间的大概走势，我们也可以推测出以后每个月的大概日照长度的大概走势。

八、第八次作业

在某化学反应中，由实验得分解物浓度与时间关系如下：

表2：分解物浓度与时间关系

时间0510152025

浓度0

时间303540455055

浓度

用最小二乘法求)

y=。

f

(t

解：利用matlab数值逼近可以得到：

即最小二乘拟合二次多项式为

-

=t

t

f

+

t

.0

.0

2453

)(2+

1992

.0

0022

其中拟合值和实际值的对比图见下

:图2：图中圆圈为已知数据的点，曲线为最小二乘拟合二次多项式

九、第九次作业

用龙贝格求积算法计算积分dx e I x ?=1

0。计算过程中数值保留6位有效数字。

1、龙贝格算法原理

龙贝格求积法的算法为

其中)(h T m 逼近I 的阶达到m 2次。

2、matlab 程序代码

function [R,quad,err,h]=romber(f,a,b,n,tol)

% f 是被积函数。

% a, b 分别为是积分的上下限。 % n ＋1是T 数表的列数。 % tol 是允许误差。 % R 是T 数表。

% quad 是所求积分值。

???

????=--==--+ ,2,1,41)

(4)2()()()(11m h T h T h T h T h T m m m m m

M=1;

h=b-a;

err=1;

J=0;

R=zeros(4,4);

R(1,1)=h*(feval('f',a)+feval('f',b))/2;

while((err>tol)&(J

J=J+1;

h=h/2;

s=0;

for p=1:M

x=a+h*(2*p-1);

s=s+feval('f',x);

end

R(J+1,1)=R(J,1)/2+h*s;

M=2*M;

for K=1:J

R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1);

end

err=abs(R(J,J)-R(J+1,K+1));%取绝对值。

end

quad=R(J+1,J+1);

再定义函数

function y=f(x);

y=exp(x);%exp(x)以e为底的指数。

命令窗口中输入

clc

clear

romber('f', 0, 1, 7, 10^(-6))

3、运行结果截图

结论：求解得最终结果为。

十、第十次作业

用经典R-K 方法计算初值问题???=-=2)0(38'y y

y ，取步长2.0=h ，计算)4.0(y 的近似值

1、经典R-K 方法设计原理

在],[1+n n x x 区间上取),(y x f 在P 个点处的数值作线性组合，即得到P 级的R-K 方法，在P 级R-K 方法中，最常用的是标准的（或经典的）四阶R-K 方法。

???

??

?

?

??

??

??++++=++=++=++==+)22(6),()2,21()

2,21(),(432113423121k k k k h y y hk y h x f k k h y h x f k k h y h x f k y x f k n n n n n

n n n n n )1.10( 2、matlab 程序代码

function fv = Funval(f,varvec,varval)

var = findsym(f);

varc = findsym(varvec); s1 = length(var); s2 = length(varc); m =floor((s1-1)/3+1); varv = zeros(1,m); if s1 ~= s2

for i=0: ((s1-1)/3)

k = findstr(varc,var(3*i+1)); index = (k-1)/3;

varv(i+1) = varval(index+1); end

fv = subs(f,var,varv); else

fv = subs(f,varvec,varval); end

function y=DELGK4_lungkuta(f, h,a,b,y0,varvec) %f:一阶常微分方程的一般表达式的右端函数 %h:积分步长

%a:自变量取值下限 %b:自变量取值上限 %y0:函数初值

%varvee ：常微分方程的变量组 format long ;

N=(b-a)/h;

y=zeros(N+1,1);

y(1)=y0;

x=a:h:b;

var=findsym(f);

for i=2:N+1

K1=Funval(f,varvec,[x(i-1) y(i-1)]);

K2=Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3=Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K2*h/2]); K4=Funval(f,varvec,[x(i-1)+h y(i-1)+h*K3]);

y(i)=y(i-1)+h*(K1+2*K2+2*K3+K4)/6;

end

format short;

3、运行结果截图

结论：最终求得)4.0(y的近似值为.|