《矩阵及数值分析》课程数值实验大作业Word格式.docx
- 文档编号:19232567
- 上传时间:2023-01-04
- 格式:DOCX
- 页数:18
- 大小:166.52KB
《矩阵及数值分析》课程数值实验大作业Word格式.docx
《《矩阵及数值分析》课程数值实验大作业Word格式.docx》由会员分享,可在线阅读,更多相关《《矩阵及数值分析》课程数值实验大作业Word格式.docx(18页珍藏版)》请在冰豆网上搜索。
formatshort
[运行结果]
4.178866707295615e-024
1/239299329230617530000000
【按第二种递推公式】
a=[11/3];
fori=2:
a=[a10/3*a(i)-a(i-1)];
数列第50项为:
2060436
【分析】
第一种算法数值稳定,计算过程舍入误差被严格控制,且按1/3的公差不断缩小。
但第二种算法数值不稳定。
另外,在第二种算法中,若将递推公式“a=[a10/3*a(i)-a(i-1)]”中的分母移动位置,改写成“a=[a10*a(i)/3-a(i-1)]”,则程序运行结果为-4966040,可以舍入误差被放大的十分严重。
二、利用迭代格式
及Aitken加速后的新迭代格式求方程
在
内的根
【未经加速的代码】
eps=1e-15;
i=1;
x0=1;
whilei<
100
x1=sqrt(10/(x0+4));
ifabs(x1-x0)<
=eps
break
end
x0=x1;
i=i+1;
方程的解[精度10^(-15)]'
disp(x1)
未经加速的迭代次数'
disp(i)
方程的解[精度10^(-15)]
1.36523001341410
未经加速的迭代次数
18
【经Aitken加速的代码】
y=sqrt(10/(x1+4));
z=sqrt(10/(y+4));
x1=z-(z-y)^2/(z-2*y+x1);
3
Aitken加速能对数列{xk}起明显的加速作用,在要求相同方程解精度的条件下,它能将迭代次数显著降低。
实际上,Aitken有时甚至能将发散的数列加速后变为收敛。
三、解线性方程组
1.分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组
迭代法计算停止的条件为:
.
2.用Gauss列主元消去法、QR方法求解如下方程组:
【1.Jacobi方法】
eps=1e-6;
A=[621-2;
250-2;
-2085;
1327];
b=[47-10]'
;
x0=zeros(4,1);
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=inv(D)*(L+U);
f=inv(D)*b;
x1=B*x0+f;
ifnorm(x1-x0)<
方程的解[精度10^(-6)]'
迭代次数'
方程的解[精度10^(-6)]
0.05204951386229
1.150********528
0.24463239101199
-0.57059207123432
迭代次数
16
【1.Gauss-Seidel方法】
B=inv(D-L)*(U);
f=inv(D-L)*b;
0.05204946628111
1.150********986
0.24463241176981
-0.57059206335719
8
【2.Gauss列主元消去法】
A=[2212;
413-1;
-4-201;
2323];
b=[1210]'
Ab=[Ab];
[N,N]=size(A);
x=zeros(N,1);
forp=1:
N-1
[max1,j]=max(abs(A(p:
N,p)));
temp=Ab(p,:
);
Ab(p,:
)=Ab(j+p-1,:
Ab(j+p-1,:
)=temp;
ifAb(p,p)==0
disp('
方程无解'
fork=p+1:
N
mult=Ab(k,p)/Ab(p,p);
Ab(k,p:
N+1)=Ab(k,p:
N+1)-mult*Ab(p,p:
N+1);
b=Ab(:
N+1);
xx=zeros(N,1);
fork=N:
-1:
1
xx(k)=b(k)/Ab(k,k);
i=(1:
k-1)'
b(i)=b(i)-xx(k)*Ab(i,k);
x=xx
x=
1.54166666666667
-2.75000000000000
.0833********
1.66666666666667
【2.QR分解法】
3
A{i}=zeros(5-i);
Q{i}=eye(4);
x=zeros(4,1);
y=zeros(4,1);
R=zeros(4);
A{1}=[2212;
413-1;
-4-201;
2323;
];
b=[1,2,1,0]'
I=eye(size(A{i}));
a=A{i}(:
1);
w=a-norm(a)*I(:
hw=I-2/(w'
*w)*(w*w'
Q{i}(i:
4,i:
4)=hw;
ifi==3
c=hw*A{i};
A{i+1}=c(2:
max(size(A{i})),2:
max(size(A{i})));
Qz=(Q{3}*Q{2}*Q{1})'
R=Q{3}*Q{2}*Q{1}*A{1};
c=Qz'
*b;
x(4)=c(4)/R(4,4);
fori=3:
x(i)=(c(i)-R(i,i+1:
4)*x(i+1:
4))/R(i,i);
x
1.54166666666666
1.66666666666666
Jacobi迭代法和Gauss-Seidel迭代法都可用来求解线性方程组。
在同等精度下,求解本道题Jacobi迭代法迭代了16次,而Gauss-Seidel仅为8次,是前者的一半。
所以可以认为,在某些情况下,Gauss-Seidel是比较好的解法,但不总如此。
Gauss消去法可能发生小主元做除数,从而导致计算结果严重偏离真实值,所以在消元过程中,每一步都按列选主元,也就是Gauss列主元消去法,它可以有效避免过于严重的舍入误差。
正交矩阵是性态最好的矩阵,它不改变矩阵的条件数,通过QR分解来计算线性方程组,也可以避免误差的放大,保证计算过程具有数值稳定性。
四、已知一组数据点
,编写一程序求解三次样条插值函数
满足
并针对下面一组具体实验数据
0.25
0.3
0.39
0.45
0.53
0.5000
0.5477
0.6245
0.6708
0.7280
求解,其中边界条件为
.
【程序代码,文件名Spline.m】
function[s]=Spline(x,y,f0,fn)
%输入实验数据x,y;
边界二阶导数f0,fn
figure
(1)
plot(x,y,'
*r'
holdon
N=max(size(x));
symsts;
fork=1:
N-1;
h(k)=x(k+1)-x(k);
g
(1)=3*(y
(2)-y
(1))/h
(1)-h
(1)/2*f0;
g(N)=3*(y(N)-y(N-1))/h(N-1)+h(N-1)/2*fn;
fork=2:
N-1
lamda(k)=h(k)/(h(k)+h(k-1));
miu(k)=h(k-1)/(h(k)+h(k-1));
g(k)=3*(miu(k)*(y(k+1)-y(k))/h(k)+...
lamda(k)*(y(k)-y(k-1))/h(k-1));
c=[1,miu(2:
N-1)];
b=2*ones(1,N);
a=[lamda(2:
N-1),1];
A=diag(c,1)+diag(b,0)+diag(a,-1);
d=c;
a=[0,a];
u
(1)=b
(1);
l(i)=a(i)/u(i-1);
u(i)=b(i)-l(i)*c(i-1);
L=eye(N)+diag(l(2:
N),-1);
U=diag(u)+diag(d,1);
z
(1)=g
(1);
z(i)=g(i)-l(i)*z(i-1);
m(N)=z(N)/u(N);
fori=N-1:
m(i)=(z(i)-c(i)*m(i+1))/u(i);
s(i)=(h(i)+2*(t-x(i)))*(t-x(i+1))^2*y(i)/(h(i))^3+...
(h(i)-2*(t-x(i+1)))*(t-x(i))^2*y(i+1)/(h(i))^3+...
(t-x(i))*(t-x(i+1))^2*m(i)/(h(i))^2+...
(t-x(i+1))*(t-x(i))^2*m(i+1)/(h(i))^2;
End
digits(8)
s=vpa(expand(s))%输出分段表达式
N-1%绘图
v=x(i):
0.005:
x(i+1);
y=subs(s(i),v);
plot(v,y);
holdon;
gridon
[命令窗口输入]
x=[0.250.30.390.450.53];
y=[0.5000,0.5477,0.6245,0.6708,0.7280];
f0=0;
fn=0;
Spline(x,y,f0,fn)
ans=
[4.6988737*t^2-.20505552*t+.35547747-6.2651650*t^3,
-2.6329843*t^2+1.9945019*t+.13552173+1.8813439*t^3,
.10638708*t^2+.92614706*t+.27440786-.45999912*t^3,
-3.4093028*t^2+2.5082075*t+.37098798e-1+2.1442156*t^3]
运行结果是一个四个元素的矩阵,各个元素依次代表四个分段区间上的三次样条曲线。
例如在第一个区间[0.250.3]上,拟合得到的三次样条曲线方程4.6988737*t^2-0.20505552*t+0.35547747-6.2651650*t^3。
由图像可知,三次样条曲线是很光滑的曲线,拟合效果很好。
五、编写程序构造区间
上的以等分结点为插值结点的Newton插值公式,假设结点数为
(包括两个端点),给定相应的函数值,插值区间和等分的份数,该程序能快速计算出相应的插值公式。
以
,
为例计算其对应的插值公式,分别取不同的
值并画出原函数的图像以及插值函数的图像,观察当
增大时的逼近效果.
【程序代码,文件名Newton.m】
functionf1=Newton(n)
symsxtf
a=-1;
b=1;
f=1./(1+25*x.^2);
y=zeros(1,n);
x=linspace(a,b,n);
h=-1:
0.01:
1;
m=n;
c(1:
m)=0.0;
yh=subs(f,h);
y=subs(f,x);
yk=y;
f1=y
(1);
y1=0;
k=1;
m-1
forj=i+1:
m
y1(j)=(y(j)-y(i))/(x(j)-x(i));
c(i)=y1(i+1);
k=k*(t-x(i));
f1=f1+c(i)*k;
y=y1;
f1=simplify(f1)
yfh=subs(f1,h);
plot(h,yfh,'
--k'
x,yk,'
x=-1:
y=1./(1+25*x.^2);
plot(x,yh,'
b'
legend('
插值曲线'
'
插值节点'
被插曲线'
Newton(4)(5个插值节点)
拟合方程为:
1-3225/754*t^2+1250/377*t^4
Newton(5)(6个插值节点)
(用digits(4),vpa(ans)处理后)
.5673+.1599e-16*t-1.731*t^2-.4441e-15*t^3+1.202*t^4+.1110e-14*t^5
当n继续增大时有:
n=7n=10
n=12n=15
本题说明在不分段拟合的条件下,节点并不是越不越好。
插值节点增多会导致插值多项式的次数增高,而高次多项式的振荡次数增多有可能使插值多项式在非节点处的误差变得很大(Runge现象)。
所以在高次插值时,为了克服这种现象,建议改为分段低次插值。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 矩阵及数值分析 矩阵 数值 分析 课程 实验 作业