矩阵实习作业.docx
《矩阵实习作业.docx》由会员分享,可在线阅读,更多相关《矩阵实习作业.docx(17页珍藏版)》请在冰豆网上搜索。
矩阵实习作业
矩阵与数值分析大作业
课程名称:
矩阵与数值分析
任课教师:
XXX(副教授)
姓名:
XXX
学号:
XXXX
班级:
XXXXX
学院:
电信与电气工程学部
交作业时间:
2011年12月21日
2011级工科硕士研究生
《矩阵与数值分析》课程数值实验题目
一、对于数列
,有如下两种生成方式给出利用上述两种递推公式生成的序列的第50项。
1、首项为
,递推公式为
;
MATLAB程序:
functionwyc11
formatlong
a=1;
fori=1:
49
1/3*a;
end
a
结果:
a=4.178866707295604e-024
2、前两项为
,递推公式为
;
MATLAB程序:
functionwyc12
formatlong
A=zeros(1,50);
A
(1)=1;
A
(2)=1/3;
fori=1:
48
A(i+2)=(10/3).*A(i+1)-A(i);
end
A(50)
结果:
a=2.060436086367324e+006
经过程序发现第二种迭代过程结果偏离正常预期结果,通过观察发现当迭代到第18项左右时,由于出现了小数减去小数的状况进而造成了误差,进而可以得出这种迭代方法不可行。
二、利用迭代格式
及Aitken加速后的新迭代格式求方程
在
内的根
首先通过普通的迭代格式求求方程
在
内的根的MATLAB程序:
functionwyc2n
c=2*10^(-6);
x=1;
x1=0;
whilec>10^(-6)
x1=sqrt(10/(x+4));
ifabs(x1-x)<10^(-6)
c=abs(x1-x);
end
x=x1;
end
x
程序结果为:
x=1.36522998713958
Aitken加速后的新迭代格式求方程
在
内的根的MATLAB程序:
functionwyc2a
formatlong
c=2*10^(-6);
x=1;
x1=0;
f=inline('sqrt(10/(x+4))','x');
whilec>10^(-6)
x1=f(f(x))-((f(f(x))-f(x))^2)/(f(f(x))-2*f(x)+x);
ifabs(x1-x)<10^(-6)
c=abs(x1-x);
end
x=x1;
end
x
程序结果为:
x=1.36523001341410
通过对比两种迭代方法可以发现Aitken加速后的新迭代格式再不出现较大误差的情况下比普通的迭代速度更快,更加的高效。
三、解线性方程组
1.分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组
迭代法计算停止的条件为:
.
Jacobi迭代法MATLAB程序:
functionwyc31j
formatlong
D=diag([6,5,8,7]);
LU=[0,-2,-1,2;
-2,0,0,2;
2,0,0,-5;
-1,-3,-2,0;];
B=inv(D)*LU;
b=[4;7;-1;0;];
f=inv(D)*b;
x=[0;0;0;0;];
x1=[0;0;0;0;];
c=2*10^(-6);
whilec>10^(-6)
x1=B*x+f;
fori=1:
3
ifabs(x1(i)-x(i))<10^(-6)
c=abs(x1(i)-x(i));
end
end
x=x1;
end
x
程序结果为:
x=[0.05205578320797
1.15094276230654
0.24463133345617
-0.57058744756179];
Gauss-Seidel迭代法MATLAB程序
functionwyc31gs
formatlong
D=diag([6,5,8,7]);
L=[0,0,0,0;
-2,0,0,0;
2,0,0,0;
-1,-3,-2,0;];
U=[0,-2,-1,2;
0,0,0,2;
0,0,0,-5;
0,0,0,0;];
b=[4;7;-1;0;];
B=inv(D-L)*U;
f=inv(D-L)*b;
x=[0;0;0;0;];
x1=[0;0;0;0;];
c=2*10^(-6);
whilec>10^(-6)
x1=B*x+f;
fori=1:
3
ifabs(x1(i)-x(i))<10^(-6)
c=abs(x1(i)-x(i));
end
end
x=x1;
end
x
程序结果为:
a=[0.05204946628111
1.150********986
0.24463241176981
-0.57059206335719
];
Gauss-Seidel迭代法对已经计算出的信息加以充分利用,可以对Jacobi迭代法进行加速提高
2.用Gauss列主元消去法、QR方法求解如下方程组:
用Gauss列主元消去法求解方程MATLAB程序:
functionwyc32g
A1=[2,2,1,2;
4,1,3,-1;
-4,-2,0,1;
2,3,2,3;];
b=[1;2;1;0];
[n,m]=size(A1);
x=zeros(n,1);
A=[A1,b];
fork=1:
n-1
forp=k+1:
n
ifabs(A(p,k))>abs(A(k,k))
forj=k:
n+1
t=A(k,j);
A(k,j)=A(p,j);
A(p,j)=t;
end
end
end
fori=k+1:
n
l=-A(i,k)/A(k,k);
forj=k+1:
n+1
A(i,j)=A(i,j)+l*A(k,j);
end
end
end
x(n)=A(n,n+1)/A(n,n);
fori=n-1:
-1:
1
s=0;
forj=i+1:
n
s=s+A(i,j)*x(j);
end
x(i)=(A(i,n+1)-s)/A(i,i);
end
x
程序结果为:
x=[1.54166666666667
-2.75000000000000
.0833********
1.66666666666667]
用QR方法求解方程MATLAB程序:
functionwyc32qr
formatlong
A=[2,2,1,2;
4,1,3,-1;
-4,-2,0,1;
2,3,2,3;];
b=[1;2;1;0];
[n,m]=size(A);
x=zeros(n,1);
Q=zeros(m,n);
R=zeros(m,n);
[Q,R]=qr(A);
b2=inv(Q)*b;
A2=[R,b2];
x(n)=A2(n,n+1)/A2(n,n);
fori=n-1:
-1:
1
s=0;
forj=i+1:
n
s=s+A2(i,j)*x(j);
end
x(i)=(A2(i,n+1)-s)/A2(i,i);
end
x
程序结果为:
x=[1.54166666666666
-2.75000000000000
.0833********
1.66666666666666]
四、已知一组数据点
,编写一程序求解三次样条插值函数
满足
并针对下面一组具体实验数据
0.25
0.3
0.39
0.45
0.53
0.5000
0.5477
0.6245
0.6708
0.7280
求解,其中边界条件为
.
MATLAB程序:
functionwyc4
x=[0.25;0.3;0.39;0.45;0.53;];
y=[0.5000;0.5477;0.6245;0.6708;0.7280;];
s21=0;
s25=0;
[n,n1]=size(x);
lamb=zeros(n,1);
mu=zeros(n,1);
h=zeros(n,1);
g=zeros(n,1);
m=zeros(n,1);
s=zeros(4,1);
symsp
fori=1:
n-1
h(i)=x(i+1)-x(i);
end
forj=2:
n-1
lamb(j)=h(j)/(h(j-1)+h(j));
mu(j)=h(j-1)/(h(j-1)+h(j));
end
lamb(n)=1;
mu
(1)=1;
fork=2:
n-1
g(k)=3*(mu(k)*(y(k+1)-y(k))/h(k)+lamb(k)*(y(k)-y(k-1))/h(k-1));
end
A=zeros(n,n);
fori=1:
n-1
A(i,i)=2;
A(i,i+1)=mu(i);
A(i+1,i)=lamb(i+1);
end
A(n,n)=2;
g
(1)=3*(y
(2)-y
(1))/h
(1)-h
(1)*s21/2;
g(n)=3*(y(n)-y(n-1))/h(n-1)+h(n-1)*s25/2;
m=inv(A)*g;
fork=1:
4
(h(k)+2.*(p-x(k))).*((p-x(k+1)).^2).*vpa(y(k)./(h(k).^3))+(h(k)-2.*(p-x(k+1))).*((p-x(k)).^2).*vpa(y(k+1)./(h(k).^3))+(p-x(k)).*((p-x(k+1)).^2).*vpa(m(k)./(h(k).^2))+(p-x(k+1)).*((p-x(k)).^2).*vpa(m(k+1)./(h(k).^2))
end
程序结果为:
当
s(p)=4000.*(-9/20+2*p)*(p-3/10)^2+4381.6*(13/20-2*p)*(p-1/4)^2+387.865164*(p-1/4)*(p-3/10)^2+369.069670*(p-3/10)*(p-1/4)^2
当
s(p)=751.303155*(-51/100+2*p)*(p-39/100)^2+856.652949*(87/100-2*p)*(p-3/10)^2+113.910392*(p-3/10)*(p-39/100)^2+98.670540*(p-39/100)*(p-3/10)^2
当
s(p)=2891.203704*(-18/25+2*p)*(p-9/20)^2+3105.555556*(24/25-2*p)*(p-39/100)^2+222.008716*(p-39/100)*(p-9/20)^2+206.2349887*(p-9/20)*(p-39/100)^2
当
s(p)=1310.15625*(-41/50+2*p)*(p-53/100)^2+1421.875*(57/50-2*p)*(p-9/20)^2+116.007181*(p-9/20)*(p-53/100)^2+109.574534*(p-53/100)*(p-9/20)^2
经过验证发现s(p)在各个区间段符合要求。
五、编写程序构造区间
上的以等分结点为插值结点的Newton插值公式,假设结点数为
(包括两个端点),给定相应的函数值,插值区间和等分的份数,该程序能快速计算出相应的插值公式。
以
,
为例计算其对应的插值公式,分别取不同的
值并画出原函数的图像以及插值函数的图像,观察当
增大时的逼近效果.
MATLAB程序:
functionwyc5
n=4;
x=sym('x','real');
A=zeros(n+1,n);
t=-1:
0.02:
1;
z=1./(1+25.*(t.^2));
f=inline('1/(1+25*(y^2))','y');
fori=1:
n+1
forj=1:
n
ifi==1
A(1,j)=2*(j-1)/(n-1)-1;
end
ifi==2
A(2,j)=f(A(1,j));
end
if(i>2)&(i ifj A(i,j)=(A(i-1,j+1)-A(i-1,j))/((i-2)*2/(n-1));
end
end
end
end
A
P3=ones(100)*A(2,1);
q=zeros(n-1,1);
fork=1:
100
x=0.02*k-1;
fori=1:
n-1
q(i)=A(i+2,1);
forj=1:
i
q(i)=q(i)*(x-A(1,i));
end
P3(k)=P3(k)+q(i);
end
end
plot(P3)
程序结果为:
当n=4时当n=6时
当n=8时当n=10时
当n=20时当n=30时
当n=50时
当n=100时
通过n取不同的值后观察发现当n的取值足够大时
区间内p(x)趋近于0,当x趋近于-1时p(x)趋近于无穷。
Newton插值公式不随插值点的增加而逼近,因此插值公式不适合对函数进行逼近。