大连理工大学
矩阵与数值分析上机作业
课程名称:矩阵与数值分析
研究生姓名:
交作业日时间:2016 年12 月20日
1.1程序:
Clear all;
n=input('请输入向量的长度n:')
for i=1:n;
v(i)=1/i;
end
Y1=norm(v,1)
Y2=norm(v,2)
Y3=norm(v,inf)
1.2结果
n=10 Y1 =2.9290
Y2 =1.2449
Y3 =1
n=100 Y1 =5.1874
Y2 =1.2787
Y3 =1
n=1000 Y1 =7.4855
Y2 =1.2822
Y3 =1
N=10000 Y1 =9.7876
Y2 =1.2825
Y3 =1
1.3 分析
一范数逐渐递增,随着n的增加,范数的增加速度减小;二范数随着n的增加,逐渐趋于定值,无群范数都是1.
2.1程序
clear all;
x(1)=-10^-15;
dx=10^-18;
L=2*10^3;
for i=1:L
y1(i)=log(1+x(i))/x(i); d=1+x(i);
if d == 1
y2(i)=1;
else
y2(i)=log(d)/(d-1);
end
x(i+1)=x(i)+dx;
end
x=x(1:length(x)-1);
plot(x,y1,'r');
hold on
plot(x,y2);
2.2 结果
2.3 分析
红色的曲线代表未考虑题中算法时的情况,如果考虑题中的算法则数值大小始终为1,这主要是由于大数加小数的原因。
第3题
3.1 程序
clear all;
A=[1 -18 144 -672 2016 -4032 5376 -4608 2304 -512];
x=1.95:0.005:2.05;
for i=1:length(x);
y1(i)=f(A,x(i));
y2(i)=(x(i)-2)^9;
end
figure(3);
plot(x,y1);
hold on;
plot(x,y2,'r');
F.m文件
function y=f(A,x) y=A(1);
for i=2:length(A); y=x*y+A(i); end;
3.2 结果
第4题
4.1 程序
clear all;
n=input('请输入向量的长度n:')
A=2*eye(n)-tril(ones(n,n),0);
for i=1:n
A(i,n)=1;
end
n=length(A);
U=A;
e=eye(n);
for i=1:n-1
[max_data,max_index]=max(abs(U(i:n,i))); e0=eye(n);
max_index=max_index+i-1;
U=e0*U;
e1=eye(n);
for j=i+1:n
e1(j,i)=-U(j,i)/U(i,i);
end
U=e1*U;
P{i}=e0;%把变换矩阵存到P中
L{i}=e1;
e=e1*e0*e;
end
for k=1:n-2
Ldot{k}=L{k};
for i=k+1:n-1
Ldot{k}=P{i}*Ldot{k}*P{i};
end
end
Ldot{n-1}=L{n-1};
LL=eye(n);
PP=eye(n);
for i=1:n-1
PP=P{i}*PP;
LL=Ldot{i}*LL;
end
b=ones(n,2);
b=e*b; %解方程
x=zeros(n,1);
x(n)=b(n)/U(n,n);
for i=n-1:-1:1
x(i)=(b(i)-U(i,:)*x)/U(i,i);
end
X=U^-1*e^-1*eye(n);%计算逆矩阵
AN=X';
result2{n-4,1}=AN;
result1{n-4,1}=x;
fprintf('%d:\n',n)
fprintf('%d ',AN);
4.2 结果
n=5
1.0625 -0.875 -0.75 -0.5 -0.0625
0.0625 1.125 -0.75 -0.5 -0.0625
0.0625 0.125 1.25 -0.5 -0.0625
0.0625 0.125 0.25 1.5 -0.0625
-0.0625 -0.125 -0.25 -0.5 0.0625
n=10
1.0625 -0.875 -0.75 -0.5 -0.0625 1.0625 -0.875 -0.75 -0.5 -0.0625 0.0625 1.125 -0.75 -0.5 -0.0625 0.0625 1.125 -0.75 -0.5 -0.0625 0.0625 0.125 1.25 -0.5 -0.0625 0.0625 0.125 1.25 -0.5 -0.0625 0.0625 0.125 0.25 1.5 -0.0625 0.0625 0.125 0.25 1.5 -0.0625 -0.0625 -0.125 -0.25 -0.5 0.0625 -0.0625 -0.125 -0.25 -0.5 0.0625 1.0625 -0.875 -0.75 -0.5 -0.0625 1.0625 -0.875 -0.75 -0.5 -0.0625 0.0625 1.125 -0.75 -0.5 -0.0625 0.0625 1.125 -0.75 -0.5 -0.0625 0.0625 0.125 1.25 -0.5 -0.0625 0.0625 0.125 1.25 -0.5 -0.0625 0.0625 0.125 0.25 1.5 -0.0625 0.0625 0.125 0.25 1.5 -0.0625 -0.0625 -0.125 -0.25 -0.5 0.0625 -0.0625 -0.125 -0.25 -0.5 0.0625
同样的方法可以算出n=20,n=30时的结果,这里就不罗列了。
第5题
5.1 程序
clear all;
n=input('请输入向量的长度n:10至20')
for i=1:n
for j=1:n
a(i,j)=1/(i+j-1);
end
end
for j=1:n
sum=0;
for k=1:j-1
sum = sum + l(j,k)^2;
end
l(j,j)=sqrt(a(j,j)-sum);
for i=j+1:n
sum=0;
for k=1:j-1
sum =sum + l(i,k)*l(j,k);
end
l(i,j)=(a(i,j)-sum)/l(j,j);
end
end
b=ones(n,1);
y=zeros(n,1);
y(n)=b(n)/l(n,n);
for i=n-1:-1:1
y(i)=(b(i)-l(i,:)*y)/l(i,i);
end
l=l';
x=zeros(n,1);
x(n)=y(n)/l(n,n);
for i=n-1:-1:1
x(i)=(y(i)-l(i,:)*x)/l(i,i);
end
fprintf('%d\t',x);
fprintf('\n');
5.2 结果
n=10 n=11 n=12 n=13 n=14 n=15 n=16 n=17 n=18 n=19 n=20
-746517.8 688 3111493.4
23
-11884558
.85
478355909
.6
497329749
.5
519377549
.9
445748685
378885969
.7
341571897
.8
102094396
994971030
.6
69820595.
08 -35465847
9.9
163411090
5
-77484610
111
-80547115
920
-82914903
907
-72327488
690
-62111481
037
-60010758
946
-1.76915E
+11
-1.68795E
+11
-15874441 97 987554340
7
-54993461
592
3.06265E+
12
3.18327E+
12
3.23484E+
12
2.86169E+
12
2.48062E+
12
2.55037E+
12
7.40163E+
12
6.93865E+
12
152439746 43 -1.17236E
+11
7.93546E+
11
-5.20373E
+13
-5.40791E
+13
-5.42806E
+13
-4.8713E+
13
-4.26882E
+13
-4.65382E
+13
-1.3123E+
14
-1.21156E
+14
-76184620 048 7.35352E+
11
-6.11155E
+12
4.7524E+1
4
4.93812E+
14
4.89563E+
14
4.46792E+
14
3.97243E+
14
4.57473E+
14
1.22E+15 1.11E+15
2.18036E+ 11 -2.70378E
+12
2.80302E+
13
-2.62E+15 -2.72E+15 -2.66E+15 -2.48E+15 -2.25E+15 -2.72E+15 -6.64E+15 -5.98E+15
-3.70513E +11 6.12295E+
12
-8.11E+13 9.24E+15 9.59E+15 9.27E+15 8.90E+15 8.30E+15 1.04E+16 2.17E+16 1.94E+16
3.69292E+ 11 -8.64269E
+12
1.51789E+
14
-2.17E+16 -2.25E+16 -2.14E+16 -2.14E+16 -2.07E+16 -2.60E+16 -4.18E+16 -3.74E+16
-1.99261E +11 7.40507E+
12
-1.83339E
+14
3.40E+16 3.53E+16 3.34E+16 3.52E+16 3.56E+16
4.28E+16 4.06E+16 3.77E+16
449059794 30 -3.52275E
+12
1.3792E+1
4
-3.55E+16 -3.68E+16 -3.47E+16 -3.89E+16 -4.12E+16 -4.38E+16 -4.89E+15 -1.01E+16 7.13565E+
11
-5.87483E
+13
2.35E+16 2.44E+16 2.37E+16 2.73E+16 2.94E+16 2.40E+16 -1.64E+16 -5.28E+15
1.08203E+
13
-8.98E+15 -9.29E+15 -1.03E+16 -9.10E+15 -7.78E+15 -4.50E+15 -1.23E+16 -1.58E+16
1.50E+15 1.55E+15 3.06E+15 -
2.94E+15 -6.68E+15 1.72E+15 1.63E+16 1.04E+16
2.17168E+
12
-7.94676E
+14
5.08E+15 7.29E+15 -1.45E+15 3.84E+16 3.85E+16
1.58892E+
14
-2.50E+15 -2.19E+15 -8.57E+15 -2.64E+16 -2.01E+16
4.80384E+
14
-3.1296E+
14
1.48E+16 -8.62E+16 -7.98E+16
2.304E+14 -9.01E+15 1.40E+17 1.19E+17
1.99E+15 -8.07E+16 -6.35E+16
1.70E+16 1.09E+16
7.5453E+1
4
第6题
6.1 程序
clear all;
A=[ 1 2 3 4;
-1 3 sqrt(2) sqrt(3);
-2 2 exp(1) pi;
-sqrt(10) 2 -3 7;
0 2 7 5/2;];
U=f61(A(:,2));
HA=f62(U,A);
f.m文件
function U=f61(x)
e1=eye(length(x),1);
U=x-sign(x(1))*sqrt(dot(x,x))*e1; U=U./sqrt(dot(U,U));
function HA=f62(U,A)
HA=A-2*U*U'*A;
6.2 结果
第7题
7.1 程序
clear all;
max=1000;
x(1,:)=[1 2 3];
for i=1:max
x(i+1,1)=(-2+x(i,2)+3*x(i,3))/5;
x(i+1,2)=(1+x(i,1)-4*x(i,3))/2;
x(i+1,3)=(10+3*x(i,1)-4*x(i,2))/15;
err(i)=sqrt(dot(x(i+1,:)-x(i,:),x(i+1,:)-x(i,:)));
if(err(i) < 10^-6)
break;
end
end
figure(7);
plot(err);
clear err;
x(1,:)=[1 2 3];
for i=1:max
x(i+1,1)=(-2+x(i,2)+3*x(i,3))/5;
x(i+1,2)=(1+x(i+1,1)-4*x(i,3))/2;
x(i+1,3)=(10+3*x(i+1,1)-4*x(i+1,2))/15;
err(i)=sqrt(dot(x(i+1,:)-x(i,:),x(i+1,:)-x(i,:)));
if(err(i) < 10^-6)
break;
end
end
hold on
plot(err,'r');
7.2 结果
误差越来越小。
第8题
8.1 程序
clear all;
max=100;
x(1)=1;
for i=1:max
x(i+1)=x(i)-(x(i)^3+2*x(i)^2+10*x(i)-100)/(3*x(i)^2+4*x(i)+10);
if(abs(x(i+1)-x(i)) < 10^-6)
break;
end
end
figure(8)
plot(x);
clear x;
x(1)=0;
x(2)=1;
for i=2:max
x(i+1)=x(i)-(x(i)^3+2*x(i)^2+10*x(i)-100)/(x(i)^3+2*x(i)^2+10*x(i)-(x (i-1)^3+2*x(i-1)^2+10*x(i-1)))*(x(i)-x(i-1));
if(abs(x(i+1)-x(i)) < 10^-6)
break;
end
end
hold on
plot(x,'r');
8.2结果
8.3 分析
由计算结果可知,弦截法的收敛速度比牛顿法的收敛速度快。第9题
9.1 程序
clear all;
f(0,4*pi);
f.m文件
function []= f( l,r )
if (r-l<10^-6)
fprintf('%g,',(r-l)/2+l);
return
end
if f9(l)*f9(r) < 0
if f9((l+r)/2+l)*f9(l)< 0
f(l,(l+r)/2+l);
end
if f9((l+r)/2+l)*f9(r)< 0
f((l+r)/2+l,r);
end
end
if f9(l)*f9(r) > 0
f(l,(l+r)/2+l);
f((l+r)/2+l,r);
end
9.2 结果
X=4.71239、8.24668、17.6715、14.1372
第10题
10.1 程序
clear all;
n=3;% 节点个数
Xj=0:1/n:1;
y=sin(pi*Xj);
for i=1:n+1
f(i,1)=y(i);
end
for j=2:n+1
for i=1:n-j+2;
dx=((j-1)/n);
f(i,j)=(f(i+1,j-1)-f(i,j-1))/dx;
end
end
for i=1:n+1
a(i)=f(1,i);
end
x=0:0.001:1;
for i=1:1/0.001+1;
y1(i)=sin(pi*x(i));
y2(i)=f10(a,Xj,x(i));
end
figure(10);
plot(x,y1,'r');
hold on;
plot(x,y2,'');
10.2结果
10.3分析
有图像可知插值函数的值已经很接近原函数的值了。第11题
11.1 程序
clear all;
n=input('请输入n:') % n代表节点
Xj=-5:1/n:5;
Yj=1./(1.+Xj.^2);
x=-5:0.01:5;
for i=1:10/0.01+1;
y1(i)=1/(1+x(i)^2);
y2(i)=f(Yj,Xj,x(i));
end
figure(11);
plot(x,y1,'r');
hold on;
plot(x,y2,,'');
f.m文件
function y=f(Yj,Xj,x)
y=0;
for i=1:length(Yj)
l=1;
for j=1:length(Xj)
if i==j
continue;
end
l=l*(x-Xj(j))/(Xj(i)-Xj(j));
end
y=y+Yj(i)*l;
end
11.2 结果
从左往右n依次取1、2、3、4
11.3 分析
随着n的不断增加,插值越来越接近真实值。
第12题
12.1 程序
clear all;
n=input('请输入n:') %n=50 100 200 500 1000
x=0:2*pi/n:2*pi;
for i=1:n+1
X=x(i);
Y(i)=exp(3*X)*cos(pi*X);
end
Y(1)=1;
for i=1:n
X=(x(i)+x(i+1))/2;
Y2(i)=exp(3*X)*cos(pi*X);
end
T=(x(n+1)-x(1))/(2*n)*(Y(1)+2*sum(Y(2:n))+Y(n+1));
S=(x(n+1)-x(1))/(6*n)*(Y(1)+4*sum(Y2)+2*sum(Y(2:n))+Y(n+1));
第13题
13.1 程序
clear all;
G2=f(1/sqrt(3))+f(-1/sqrt(3));
G3=0.555555556*f(-0.7745966692)+0.555555556*f(+0.7745966692)+0.888888 88889*f(0);
G5=0.2369268851*f(-0.9061798459)+0.2369268851*f(+0.9061798459)+...
0.4786286705*f(-0.5384693101)+0.4786286705*f(+0.5384693101)+...0.5688 888889*f(0);
f.m文件1
function y=f(x)
y=x^2/sqrt(1-x^2);
f.m文件2
function y=f(x)
x=pi/2*(x+1)/2;
y=sin(x)/x;
13.2
第14题
14.1 程序
clear all;
U0=2;%初值
step=0.2;%step表示步长
U1(1)=U0;
Len=1;
for i=1:floor(Len/step)
U1(i+1)=U1(i)+step*f14(U1(i),(i-1)*step);
end
% 改进的Euler法
U2(1)=U0;
for i=1:floor(Len/step)
U2(i+1)=U2(i)+step/2*(f14(U2(i),(i-1)*step)+f14(U1(i+1),(i+1-1)*step) );
end
% Runge_KUtta
U3(1)=U0;
for i=1:floor(Len/step)
t=(i-1)*step;
k1=f14(U3(i),t);
k2=f14(U3(i)+step/2*k1,t+step/2);
k3=f14(U3(i)+step/2*k2,t+step/2);
k4=f14(U3(i)+step*k3,t+step);
U3(i+1)=U3(i)+step/6*(k1+2*k2+2*k3+k4); end
for i=1:floor(Len/step)+1
t(i)=(i-1)*step;
end
figure(14)
plot(t,U1,'');
plot(t,U2,'g');
plot(t,U3,'r');
f.m文件
function y=f(x,t)
y=1/((t+1)^2)*(t*x-x^2);
14.2 结果
数值分析上机实验报告 选题:曲线拟合的最小二乘法 指导老师: 专业: 学号: 姓名:
课题八曲线拟合的最小二乘法 一、问题提出 从随机的数据中找出其规律性,给出其近似表达式的问题,在生产实践和科学实验中大量存在,通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中,根据统计数据的含碳量与时间关系,试求含碳量y 与时间t 的拟合曲线。 二、要求 1、用最小二乘法进行曲线拟合; 2、近似解析表达式为()33221t a t a t a t ++=?; 3、打印出拟合函数()t ?,并打印出()j t ?与()j t y 的误差,12,,2,1 =j ; 4、另外选取一个近似表达式,尝试拟合效果的比较; 5、*绘制出曲线拟合图*。 三、目的和意义 1、掌握曲线拟合的最小二乘法; 2、最小二乘法亦可用于解超定线代数方程组; 3、探索拟合函数的选择与拟合精度间的关系。 四、计算公式 对于给定的测量数据(x i ,f i )(i=1,2,…,n ),设函数分布为 ∑==m j j j x a x y 0)()(? 特别的,取)(x j ?为多项式 j j x x =)(? (j=0, 1,…,m )
则根据最小二乘法原理,可以构造泛函 ∑∑==-=n i m j i j j i m x a f a a a H 1 10))((),,,(? 令 0=??k a H (k=0, 1,…,m ) 则可以得到法方程 ???? ??????? ?=????????????????????????),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ????????????????????? 求该解方程组,则可以得到解m a a a ,,,10 ,因此可得到数据的最小二乘解 ∑=≈m j j j x a x f 0)()(? 曲线拟合:实际工作中,变量间未必都有线性关系,如服药后血药浓度与时间的关系;疾病疗效与疗程长短的关系;毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据,并用拟合的曲线方程分析两变量间的关系。 五、结构程序设计 在程序结构方面主要是按照顺序结构进行设计,在进行曲线的拟合时,为了进行比较,在程序设计中,直接调用了最小二乘法的拟合函数polyfit ,并且依次调用了plot 、figure 、hold on 函数进行图象的绘制,最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差,进行拟合效果的比较。
工科数学分析上机作业 说明:以下两道题均是使用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
昆明理工大学工科研究生《数值分析》上机实验 学院:材料科学与工程学院 专业:材料物理与化学 学号:2011230024 姓名: 郑录 任课教师:胡杰
P277-E1 1.已知矩阵A= 10787 7565 86109 75910 ?? ?? ?? ?? ?? ??,B= 23456 44567 03678 00289 00010 ?? ?? ?? ?? ?? ?? ?? ?? ,错误!未找到引用源。 = 11/21/31/41/51/6 1/21/31/41/51/61/7 1/31/41/51/61/71/8 1/41/51/61/71/81/9 1/51/61/71/81/91/10 1/61/71/81/91/101/11?????????????????? (1)用MA TLAB函数“eig”求矩阵全部特征值。 (2)用基本QR算法求全部特征值(可用MA TLAB函数“qr”实现矩阵的QR分解)。解:MA TLAB程序如下: 求矩阵A的特征值: clear; A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; E=eig(A) 输出结果: 求矩阵B的特征值: clear; B=[2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0]; E=eig(B) 输出结果:
求矩阵错误!未找到引用源。的特征值: clear; 错误!未找到引用源。=[1 1/2 1/3 1/4 1/5 1/6; 1/2 1/3 1/4 1/5 1/6 1/7; 1/3 1/4 1/5 1/6 1/7 1/8; 1/4 1/5 1/6 1/7 1/8 1/9;1/5 1/6 1/7 1/8 1/9 1/10; 1/6 1/7 1/8 1/9 1/10 1/11]; E=eig(错误!未找到引用源。) 输出结果: (2)A= 10 7877565861097 5 9 10 第一步:A0=hess(A);[Q0,R0]=qr(A0);A1=R0*Q0 返回得到: 第二部:[Q1,R1]=qr(A1);A2=R1*Q1
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);
2015.1.9 上机作业题报告 USER
1.Chapter1 1.1题目 设S N = 1 j 2?1 N j =2 ,其精确值为 )1 1123(21+--N N 。 (1)编制按从大到小的顺序1 1 1311212 22-+??+-+-=N S N ,计算S N 的通用程序。 (2)编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ,计算S N 的通用程序。 (3)按两种顺序分别计算64210,10,10S S S ,并指出有效位数。(编制程序时用单精度) (4)通过本次上机题,你明白了什么? 1.2程序 1.3运行结果
1.4结果分析 按从大到小的顺序,有效位数分别为:6,4,3。 按从小到大的顺序,有效位数分别为:5,6,6。 可以看出,不同的算法造成的误差限是不同的,好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时,误差限随着N 的增大而增大,可见在累加的过程中,误差在放大,造成结果的误差较大。因此,采取从小到大的顺序累加得到的结果更加精确。 2.Chapter2 2.1题目 (1)给定初值0x 及容许误差ε,编制牛顿法解方程f(x)=0的通用程序。 (2)给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321= *=*-=*x x x ○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时,Newton 迭代序列收敛于根x2*。试确定尽可能大的δ。 ○2试取若干初始值,观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 (3)通过本上机题,你明白了什么? 2.2程序
东南大学数值分析上机作业 汇总 -标准化文件发布号:(9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII
数值分析上机报告 院系: 学号: 姓名:
目录 作业1、舍入误差与有效数 (1) 1、函数文件cxdd.m (1) 2、函数文件cddx.m (1) 3、两种方法有效位数对比 (1) 4、心得 (2) 作业2、Newton迭代法 (2) 1、通用程序函数文件 (3) 2、局部收敛性 (4) (1)最大δ值文件 (4) (2)验证局部收敛性 (4) 3、心得 (6) 作业3、列主元素Gauss消去法 (7) 1、列主元Gauss消去法的通用程序 (7) 2、解题中线性方程组 (7) 3、心得 (9) 作业4、三次样条插值函数 (10) 1、第一型三次样条插值函数通用程序: (10) 2、数据输入及计算结果 (12)
作业1、舍入误差与有效数 设∑ =-=N j N j S 2 2 11 ,其精确值为?? ? ??---1112321N N . (1)编制按从小到大的顺序1 1 131121222-? ??+-+-=N S N ,计算N S 的通用程序; (2)编制按从大到小的顺序()1 21 11111222-???+--+-=N N S N ,计算N S 的通用程序; (3)按两种顺序分别计算642101010,,S S S ,并指出有效位数; (4)通过本上机你明白了什么? 程序: 1、函数文件cxdd.m function S=cxdd(N) S=0; i=2.0; while (i<=N) S=S+1.0/(i*i-1); i=i+1; end script 运行结果(省略>>): S=cxdd(80) S= 0.737577 2、函数文件cddx.m function S=cddx (N) S=0; for i=N:-1:2 S=S+1/(i*i-1); end script 运行结果(省略>>): S=cddx(80) S= 0.737577 3、两种方法有效位数对比
函数定义: % 目标函数 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;
¥ 数值分析思考题1 1、讨论绝对误差(限)、相对误差(限)与有效数字之间的关系。 2、相对误差在什么情况下可以用下式代替 3、查阅何谓问题的“病态性”,并区分与“数值稳定性”的不同点。 4、取 ,计算 ,下列方法中哪种最好为什么(1)(3 3-,(2)(2 7-,(3) ()3 1 3+ ,(4) ()6 1 1 ,(5)99- , 数值实验 数值实验综述:线性代数方程组的解法是一切科学计算的基础与核心问题。求解方法大致可分为直接法和迭代法两大类。直接法——指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法,因此也称为精确法。当系数矩阵是方的、稠密的、无任何特殊结构的中小规模线性方程组时,Gauss消去法是目前最基本和常用的方法。如若系数矩阵具有某种特殊形式,则为了尽可能地减少计算量与存储量,需采用其他专门的方法来求解。 Gauss消去等同于矩阵的三角分解,但它存在潜在的不稳定性,故需要选主元素。对正定对称矩阵,采用平方根方法无需选主元。方程组的性态与方程组的条件数有关,对于病态的方程组必须采用特殊的方法进行求解。 数值计算方法上机题目1 1、实验1. 病态问题 实验目的: 算法有“优”与“劣”之分,问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感,反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中,如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决,当然一般要付出一些代价(如耗用更多的机器时间、占用更多的存储空间等)。 $ r e x x e x x ** * ** - == 141 . ≈)61