大连理工大学矩阵与数值分析上机作业
————————————————————————————————作者：————————————————————————————————日期：
矩阵与数值分析上机作业
学校：大连理工大学
学院：
班级：
姓名：
学号：
授课老师：注：编程语言Matlab
程序：
Norm.m函数
function s=Norm(x,m)
%求向量x的范数
%m取1,2,inf分别表示1,2，无穷范数n=length(x);
s=0;
switch m
case 1 %1-范数
for i=1:n
s=s+abs(x(i));
end
case 2 %2-范数
for i=1:n
s=s+x(i)^2;
end
s=sqrt(s);
case inf %无穷-范数
s=max(abs(x));
end
计算向量x，y的范数
Test1.m
clear all;
clc;
n1=10;n2=100;n3=1000;
x1=1./[1:n1]';x2=1./[1:n2]';x3=1./[1:n3]'; y1=[1:n1]';y2=[1:n2]';y3=[1:n3]';
disp('n=10时');
disp('x的1-范数:');disp(Norm(x1,1));
disp('x的2-范数:');disp(Norm(x1,2));
disp('x的无穷-范数:');disp(Norm(x1,inf));
disp('y的1-范数:');disp(Norm(y1,1));
disp('y的2-范数:');disp(Norm(y1,2));
disp('y的无穷-范数:');disp(Norm(y1,inf));
disp('n=100时');
disp('x的1-范数:');disp(Norm(x2,1));
disp('x的2-范数:');disp(Norm(x2,2));
disp('x的无穷-范数:');disp(Norm(x2,inf));
disp('y的1-范数:');disp(Norm(y2,1));
disp('y的2-范数:');disp(Norm(y2,2));
disp('y的无穷-范数:');disp(Norm(y2,inf));
disp('n=1000时');
disp('x的1-范数:');disp(Norm(x3,1));
disp('x的2-范数:');disp(Norm(x3,2));
disp('x的无穷-范数:');disp(Norm(x3,inf));
disp('y的1-范数:');disp(Norm(y3,1));
disp('y的2-范数:');disp(Norm(y3,2));
disp('y的无穷-范数:');disp(Norm(y3,inf));
运行结果：
n=10时
x的1-范数:2.9290；x的2-范数:1.2449； x的无穷-范数:1
y的1-范数:55； y的2-范数:19.6214； y的无穷-范数:10
n=100时
x的1-范数:5.1874；x的2-范数: 1.2787； x的无穷-范数:1
y的1-范数:5050； y的2-范数:581.6786； y的无穷-范数:100
n=1000时
x的1-范数:7.4855； x的2-范数:1.2822； x的无穷-范数:1
y的1-范数: 500500； y的2-范数:1.8271e+004；y的无穷-范数:1000
程序
Test2.m
clear all;
clc;
n=100;%区间
h=2*10^(-15)/n;%步长
x=-10^(-15):h:10^(-15);
%第一种原函数
f1=zeros(1,n+1);
for k=1:n+1
if x(k)~=0
f1(k)=log(1+x(k))/x(k);
else
f1(k)=1;
end
end
subplot(2,1,1);
plot(x,f1,'-r');
axis([-10^(-15),10^(-15),-1,2]); legend('原图');
%第二种算法
f2=zeros(1,n+1);
for k=1:n+1
d=1+x(k);
if(d~=1)
f2(k)=log(d)/(d-1);
else
f2(k)=1;
end
end
subplot(2,1,2);
plot(x,f2,'-r');
axis([-10^(-15),10^(-15),-1,2]); legend('第二种算法');
运行结果：
显然第二种算法结果不准确，是因为计算机中的舍入误差造成的，当]10,10[1515--∈x 时，
x d +=1，计算机进行舍入造成d 恒等于1，结果函数值恒为1。
程序：
秦九韶算法:QinJS.m
function y=QinJS(a,x) %y 输出函数值
%a 多项式系数，由高次到零次 %x 给定点 n=length(a); s=a(1); for i=2:n
s=s*x+a(i);
end
y=s;
计算p(x):test3.m
clear all;
clc;
x=1.6:0.2:2.4;%x=2的邻域
disp('x=2的邻域：');x
a=[1 -18 144 -672 2016 -4032 5376 -4608 2304 -512]; p=zeros(1,5);
for i=1:5
p(i)=QinJS(a,x(i));
end
disp('相应多项式p值：');p
xk=1.95:0.01:20.5;
nk=length(xk);
pk=zeros(1,nk);
k=1;
for k=1:nk
pk(k)=QinJS(a,xk(k));
end
plot(xk,pk,'-r');
xlabel('x');ylabel('p(x)');
数值分析上机题 第一章 17.（上机题）舍入误差与有效数 设∑=-= N j N j S 2 2 11 ，其精确值为）111-23（21+-N N 。 （1）编制按从大到小的顺序1 -1 ···1-311-21222N S N +++=，计算N S 的通用 程序； （2）编制按从小到大的顺序1 21 ···1)1(111 222-++--+ -=N N S N ，计算N S 的通用程序； （3）按两种顺序分别计算210S ,410S ,610S ,并指出有效位数（编制程序时用单精度）； （4）通过本上机题，你明白了什么？ 解： 程序： （1）从大到小的顺序计算1 -1 ···1-311-21222N S N +++= ： function sn1=fromlarge(n) %从大到小计算sn1 format long ; sn1=single(0); for m=2:1:n sn1=sn1+1/(m^2-1); end end （2）从小到大计算1 21 ···1)1(111 2 22 -++--+-= N N S N function sn2=fromsmall(n) %从小到大计算sn2 format long ; sn2=single(0); for m=n:-1:2 sn2=sn2+1/(m^2-1); end end （3） 总的编程程序为： function p203()
clear all format long; n=input('please enter a number as the n:') sn=1/2*(3/2-1/n-1/(n+1));%精确值为sn fprintf('精确值为%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:n sn1=sn1+1/(m^2-1); end end function sn2=fromsmall(n) %从小到大计算sn2 format long; sn2=single(0); for m=n:-1:2 sn2=sn2+1/(m^2-1); end end end 运行结果：
工科数学分析上机作业 说明：以下两道题均是使用Matlab 语言，且在Matlab 7.0中运行通过。 1.（两个重要极限）计算下列函数的函数值并画出图形，观察两个重要极限值。 (1)y=f(x)=; (2)y=f(x)=. sin x x (1+x)1x 解：（1）求解过程如下： >> syms x >> y=limit(sin(x)/x) y = 1 >> ezplot(sin(x)/x,[-10*pi,10*pi]) >> ezplot(sin(x)/x,[-1*pi,1*pi]) 其图形如下：
（2）求解过程如下：>> syms x >> y=(1+x)^(1/x)
y = (1+x)^(1/x) >> y=limit((1+x)^(1/x)) y = exp(1) >> ezplot((1+x)^(1/x),[-1000,1000]) >> ezplot((1+x)^(1/x),[-10,10]) >> ezplot((1+x)^(1/x),[-1,1]) 其图像如下：
分析如下：（1）当x 取值为[-30,30]时，由该题的第一个图像可以看到，函数值在不断震荡，一会为正数，一会为负数。
而当x 取值为[-3，3]时，函数值始终大于0。当x 趋近于0时，由该题的第二个图像可以得到函数值为1。 另外，该结论也可以由夹逼法则证明，结果不变，当x 趋近于0时，函数值仍为1。 （2）由该题的三个图像可以知道，该函数在定义域内为单调递减函数。且由该题的第一和二个图像知道，当x 在 [0,10]区间内，函数递减趋势非常迅速。由该题的第三个图像知道，当x 趋于0 时，函数值为自然对数的底数 e ,即约为2.71828. 3.计算f(x)=, 12+1√2π ∫x 0e ?t 2/2dt 1?x ?3的函数值{f （0.1k ）;k=1,2,…,30}.计算结果取7位有效数字。 解：计算过程为： >> f1= @(t) exp(-(t).^2/2) f1 = @(t) exp(-(t).^2/2) >> for i=1:30
实验一误差分析 实验1.1（病态问题） 实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出：考虑一个高次的代数多项式 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 其中ε（1.1）和（1.221,,,a a 的输出b ”和“poly ε。 （1（2 （3）写成展 关于α solve 来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考Matlab 的帮助。 实验过程： 程序： a=poly(1:20); rr=roots(a); forn=2:21 n form=1:9 ess=10^(-6-m);
ve=zeros(1,21); ve(n)=ess; r=roots(a+ve); -6-m s=max(abs(r-rr)) end end 利用符号函数：（思考题一）a=poly(1:20); y=poly2sym(a); rr=solve(y) n
很容易的得出对一个多次的代数多项式的其中某一项进行很小的扰动，对其多项式的根会有一定的扰动的，所以对于这类病态问题可以借助于MATLAB来进行问题的分析。 学号：06450210 姓名：万轩 实验二插值法
2016年理工大学优化方法上机大作业学院： 专业： 班级： 学号： ： 上机大作业1： 1.最速下降法：
function f = fun(x) f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; end function g = grad(x) g = zeros(2,1); g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2); end
function x_star = steepest(x0,eps) gk = grad(x0); res = norm(gk); k = 0; while res > eps && k<=1000 dk = -gk; ak =1; f0 = fun(x0); f1 = fun(x0+ak*dk); slope = dot(gk,dk); while f1 > f0 + 0.1*ak*slope ak = ak/4; xk = x0 + ak*dk; f1 = fun(xk); end k = k+1; x0 = xk; gk = grad(xk); res = norm(gk); fprintf('--The %d-th iter, the residual is %f\n',k,res); end x_star = xk; end >> clear
>> x0=[0,0]'; >> eps=1e-4; >> x=steepest(x0,eps)
2.牛顿法： function f = fun(x) f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; end function g = grad2(x) g = zeros(2,2);
数值分析上机实验报告
《数值分析》上机实验报告 1.用Newton 法求方程 X 7-X 4+14=0 在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。 1.1 理论依据： 设函数在有限区间[a ，b]上二阶导数存在，且满足条件 {}α?上的惟一解在区间平方收敛于方程所生的迭代序列 迭代过程由则对任意初始近似值达到的一个中使是其中上不变号 在区间],[0)(3,2,1,0,) (') ()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20 )()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f a b c f x f b a x f b f x f k k k k k k ==- ==∈≤-≠>+ 令 )9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3 2 2 5 333647>?''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f 故以1.9为起点 ?? ?? ? ='- =+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代，逼近x 的真实根。当前后两个的差<=ε时，就认为求出了近似的根。本程序用Newton 法求代数方程(最高次数不大于10)在（a,b ）区间的根。
1.2 C语言程序原代码： #include
函数定义： % 目标函数 function f = fun(x) fm=0; for i=1:499 fmi = (1-x(i))^2 + 100*(x(i+1)-x(i)^2)^2; fm=fm+fmi; end f =fm; end % 梯度 function g = grad(x) g = zeros(500,1); g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); for i=2:499 g(i)=2*(x(i)-1)+400*x(i)*(x(i)^2-x(i+1))+200*(x(i)-x(i-1)^2); end g(500) = 200*(x(500)-x(499)^2); end % 二阶梯度
function g = grad2(x) g = zeros(500,500); g(1,1)=2+400*(3*x(1)^2-x(2)); g(1,2)=-400*x(1); for i=3:500 g(1,i)=0; end for i=1:498 g(500,i)=0; end g(500,499)=-400*x(499); g(500,500)=200; for i=2:499 for j=1:500 if j==i-1 g(i,j)= -400*x(i-1); elseif j==i g(i,j)= 2+400*(3*x(i)^2-x(i+1))+200; elseif j==i+1 g(i,j)= -400*x(i); else g(i,j)=0; end end end end 1.最速下降法 function x_star = steepest(x0,eps) gk = grad(x0); res = norm(gk); k = 0; while res > eps && k<=50000 dk = -gk;