版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
矩阵与数值分析实验报告
任课教师:________张宏伟_________________________
教学班号:__________02___________________________
院系:电信(计算机应用技术)
学生姓名:
学生学号:
1.方程,f(c)=2ir3-5rr2-19a:+42=0在x=3.0附近有根,试写出其三种不同的
等价形式以构成两种不同的迭代格式,再用这两种迭代求根,并绘制误差下
降曲线,观测这两种迭代是否收敛及收敛的快慢。
解答:
代码如下:
clear;
symsxty;
x=3;%初始迭代点
t=0;%中间变量
y=o;%绘制下降曲线变化量
k=0;/迭代计数变量
epx=1;%变量计算差
E=le-20;%精度
f1=(2*x~3-5*x-2+42)/19;%迭代1
f2=(2*x-3-19*x+42)”(l/2);%迭代2
f3=(5*x*2+19*x-42)*(1/3);%迭代3
whi1e(epx>E&k<1000)%循环迭代
k=k+l;
y(k)=x;
t=f1;%标记1
epx=abs(t-x);
x=t;
f1=(2*x*3-5*x,2+42)/19;%标记2
end;
Plot(y);
title('迭代1误差下降曲线');
迭代公式1收敛结果:x=2
迭代公式1误差变化曲线
x10迭代2误差下降曲线
3.5
2.5
2
1.5
0.5
10
迭代公式2误差变化曲线
迭代公式3收敛结果:x=6.8750
迭代公式3误差变化曲线
结果分析:迭代公式1:fl=(2*x*3-5*x*2+42)/19和迭代公式3:f3=(5*x-2
+19*x-42)71/3)计算结果收敛,其中迭代公式1收敛速度快于迭代公式3o
迭代公式2:f2=(2*x*3-19*x+42)71/2)计算结果不收敛。
2.用复化梯形公式、复化辛普森公式、龙贝格公式求下列定积分,规定绝对误
差为,=().5x10-8并将计算结果与精确解进行比较:
4232
e=-£e.2n6=
3
解答:
(1)复化梯形公式代码:
clear;
symsxsum1sum2albl;
fun=(2/3)*(x-3)*exp(x-2);
suml=0;%积分结果变量1
sum2=0;%积分结果变量2
n=1;%迭代计数变量
epx=l阶段误差
E=0.5e-8;%精度
a=l;%积分下限
b=2;%积分上限
while(epx>E)
h=(b—a)/n;
fori=0:(nT)%for循环求解一次N点积分结果
a1=subs(fun,x,(a+i*h));
bl=subs(fun,x,(a+(i+1)*h));
t=((al+bl)*h)/2;
sum2=sum2+t;
end;
epx=abs(suml-sum2);%计算阶段误差
suml=sum2;为使用suml暂存上次的计算结果
sum2=0;
n=n*2;
end;
disp('积分结果为:,),
vpa(sum1)
积分结果为:
ans=
54.4316975
真值约为:54.598
分析:由结果可知计算结果具有8位有效数字,已满足精度规定。使用数据点数
为N=8192o
复化辛普森公式代码:
clear;
symsxsumlsum2a1b1;
fun=(2/3)*(x-3)*exp(x'2);
suml=0;%积分结果变量1
sum2=0;%积分结果变量2
n=l;%迭代计数变量
epx=1;%变量计算差
E=0.5e-8;%精度
a=l;%积分下限
b=2;%积分上限
whi1e(epx>E)
h=(b—a)/n;
fori=0:(n-1)
a1=subs(fun,x,(a+i*h));
bl=subs(fun,x,(a+(i+1)*h));
ab=subs(fun,x,((a+i*h)+(a+(i+l)*h))/2);
t=((al+4*ab+bl)*h)/6;
sum2=sum2+t;
end;
epx=abs(sum1-sum2);
suml=sum2;
sum2=0;
n=n*2;
end;
disp('积分结果为:’),
vpa(suml)
积分结果为:
ans=
54.59858433
真值约为:54.598
分析:由结果可知计算结果具有11位有效数字,已满足精度规定。使用数据点
数为N=1024,对比复化梯形公式可以复化辛普森公式计算精度更高。
龙贝格公式代码:
clear;
symsxa1bly;
fun=(2/3)*(x-3)*exp(x-2);
suml=O;96积分结果变量1
sum2=0;%积分结果变量2
epx=l;%变量计算差
E=0.5e-8;蟒青度
a=1;%积分下限
b=2;%积分上限
k=0;%迭代计数变量1
m=0;%迭代计数变量2
while(epx>E)
k=k+1;
m=m+1;
h=(b—a)/(2-k);
fori=0:⑵k—l)
al=subs(fun,x,(a+i*h));
b1=subs(fun,x,(a+(i+1)*h));
t=((al+bl)*h)/2;
sum2=sum2+t;
end;
sum1=sum2;
y(m)=sum2;
sum2=0;
fori=0:k-2;
m=m+1;
y(m)=((4"(k—1))*y(m—1)—y(m—k))/(4^(k-l)—1);
end;
if(m>l)
epx=abs(sym(y(m)-y(m-1)));
end;
end;
disp(,积分结果为:'),
Vpa(y(m))
积分结果为:
ans=
54.6804633
真值约为:54.598
分析:由结果可知计算结果具有9位有效数字,已满足精度规定。
(2)积分结果为:
ans=
1.96865767
真值约为:1.
分析:由结果可知计算结果具有8位有效数字,已满足精度规定。使用数据点数为
N=65536O
积分结果为:
ans=
1.48476064
真值约为:1.
分析:由结果可知计算结果具有10位有效数字,已满足精度规定。使用数据点数
为N=512,可知复化辛普森公式精度明显高于复化梯形公式。
积分结果为:
ans=
1.
真值约为:1.
分析:由结果可知计算结果具有8位有效数字,已满足精度规定。
3.使用带选主元的分解法求解线性方程组AT=h,其中4=S’lnx",3=",
^ilnxi,bn=n.当,>2时,bn=(*-1)/("1).对于"=37”的情况分别求解。
精确解为工=(11---IIf对得到的结果与精确解的差异进行解释。
解答:
程序代码:
functionx=luso1ve(A,b)
[n,n]=size(A);
x=zeros(n,1);
y=zeros(n,1);
temprow=zeros(n,1);
tempconstant=0;
Pvector=zeros(n,1);
forcol=l:n-1
[max_element,index]=max(abs(A(co1:n,col)));
temprow=A(col,:);
A(col,:)=A(index+co1-1,:);
A(index+col-1,:)=temprow;
tempconstant=b(co1);
b(col)二b(index+co1—1);
b(index+col-1)=tempconstant;
ifA(co1,col)=0
disp('Aissingular.nouniquesolution');
return
end
forrow=col+1:n
mult=A(row,col)/A(col,col);
A(row,col)=mult;
A(row,col+1:n)=A(row,col+1:n)—mul
t*A(col,col+1:n);
end
end
y(1)=b(l);
fork=2:n
y(k)=b(k)—A(k,1:k-1)*y(l:k-1);
end
x(n)=y(n)/A(n,n);
fork=n-1:-l:1
x(k)=(y(k)-A(k,k+1:n)*x(k+l:n))/A(k,k);
end
实验结果:
T
n=3时,解得x=[1,1,1];
T
n=7时,解得x=[1,1,1,1,1>1,1];
n=11时,解得x=[1.0001,0.9998,1.0002,0.9999,1,1,1,1,1,
1,I;
实验结果分析:
高斯列主元消去法,避免了小主元做除数,在Gauss消去法中增长选主
元的过程,在第K步消元时,一方面在第K列主对角元以下元素中挑选绝对值
最大的数,并通过初等行变换,使得该数位于主对角上,然后继续消元。当矩
阵维数较小时,精度较高,但随着矩阵维数增大,精度下降。
4.用4阶Runge-kutta法求解微分方程
u=e~2'-2u,u(0)=—,u(t)=-e-%+te2'
1010
(1)令也=0],使用上述程序执行20步,然后令人=005,使用上述程序执行4
0步
(2)比较两个近似解与精确解
(3)当〃减半时,(1)中的最终全局误差是否和预期相符?
(4)在同一坐标系上画出两个近似解与精确解.(提醒输出矩阵R包含近似解的
x和V坐标,用命令pl。t(R(:,1),R(:,2))画出相应图形.)
解答:
程序代码:
functionR=rk4(f,a,b,ya,n)
symstu
h=(b—a)/n;
T=zeros(1,n+1);
Y=zeros(l,n+1);
T=a:h:b;
Y(1)=ya;
fork=1:n
K1=h*subs(f,{t,u},{T(k),Y(k)});
K2=h*subs(f,{t,u},{T(k)+h/2,Y(k)+Kl/2});
K3=h*subs(f,{t,u},{T(k)+h/2,Y(k)+K2/2});
K4=h*subs(f,{t,u},{T(k)+h,Y(k)+K3));
Y(k+1)=Y(k)+(K1+2*K2+2*K3+K4)/6;
end
R[T'Y'];
h=0.1时,得到的解:
t近似解精确解
00.000.00
0.000.2800.560
0.000.470.69
0.000.2230.261
0.000.870.61
0.000.090.87
0.000.270.54
0.000.0970.329
0.000.5480.519
0.000.9920.159
1.000.8250.027
1.000.8460.480
1.000.1770.124
1.000.410.07
1.000.090.83
1.000.490.58
1.000.720.22
1.000.820.59
1.000.0500.056
1.000.010.33
2.000.640.34
h=0.05时得到的解:
t近似解精确解
00.000.00
0.000.8840.539
0.000.1210.560
0.0000.3970.043
0.000.430.69
0.000.870.42
0.000.2190.261
0.000.530.13
0.000.820.61
0.000.570.33
0.000.860.87
0.000.2350.275
0.000.210.54
0.000.600.51
0.000.2780.329
0.000.0230.617
0.000.6380.519
0.000.5120.010
0.000.7010.159
0.000.4030.377
1.000.9070.027
1.000.7040.093
1.000.3140.480
1.0000.100.50
1.000.1850.124
1.000.100.26
1.000.120.07
1.000.010.64
1.000.370.83
1.000.030.43
1.000.300.58
1.000.690.37
1.000.530.22
1.000.660.17
1.000.700.59
1.000.770.29
1.000.0590.056
1.000.780.16
1.000.660.33
1.000.0400.040
2.000.420.34
当h减半时,(1)中的全局误差缩小,和最终的预期相符。
近似解与精确解图形:
5.设
2
-1
Tn=
-1
2
为”阶的三对角方阵,4,是一个〃阶的对称正定矩阵
Tn+21n—1n
—InTn+27n
~In
—InTn+2In
其中心为加介单位矩阵。设工=(1.1,1,…,1)7为线性方程组=/>的真解,右边的
向量/曲这个真解给出。
(1)用Cholesky分解法求解该方程。
⑵用Jacobi迭代法和Gauss-Seidel迭代法求解该方程组,误差设为
E=IN
其中”取值为4,5,6.
解答:
(l)Cholesky分解程序代码(Ma11ab)
只列出n等于4时的代码,n等于5和6时代码类似。
c1ear
clc
%定义T(n)
T4=[2-100;-12-10;0-12-1;00-12];
%定义单位矩阵
14=eye(4);
/定义0矩阵
Z4=zeros(4,4);
/定义A(n)
%定义A(n)中的对角元素
d4=T4+2*14;
A4=[d4-14Z4Z4;-I4d4-14Z4;Z4-I4d4-I4;Z4Z4-
14d4];
%定义x
x4=ones(16,1);
%求右边向量b
b4=A4*x4;
%(1)用Cho1esky分解法求解方程组
14=zeros(16,16);
%n为4时求解L
n=16;
forj=l:n
sum=0;
fork=l:j-1
sum=sum+l4(j,k)14(j,k);
end
14(j,j)=sqrt(A4(j,j)-sum);
fori=j+l:n
sum=0;
fork=l:j-1
sum=sum+14(i,k)*14(j,k);
end
14(i,j)=(A4(i,j)—sum)/14(j,j)
end
end
%n等于4时求解方程的解LL'x=b
%Ly=bL'x=y
%计算y
y4=zeros(16,1);
y4(l)=b4(1)/14(1,1);
fori=1:16
sum=0;
fork=1:i—1
sum=sum+l4(i,k)*y4(k);
end
y4(i)=(b4(i)-sum)/l4(i,i);
end
%计算x_4
x_4=zeros(16,1);
x_4(l6)=y4(l6)/14(16,16);
fori=15:-l:l
sum=0;
fork=i+l:16
sum=sum+14(k,i)*x_4(k);
end
x_4(i)=(y4(i)-sum)/14(i,i);
end
运营结果:(n等于4,5,6时)x的转置
n=4时:[1111111111111111]
n=5时:[1111111111111111111111
1111
n=6时:[111111111111111111111
111111111111111]
(2)Jacobi迭代法代码(只列出n等于4时的代码,n等于5和6时类似)
%取初始向量为0
jx40=zeros(16,1);
jx41=zeros(16,1);
dx4=zeros(16,1);
while1
fori=l:16
sum=0;
forj=l:16
sum=sum+A4(i,j)*jx40(j);
end
dx4(i)=(b4(i)—sum)/A4(i,i);
jx41(i)=jx40(i)+dx4(i);
end
ifnorm(jx41-jx40)>le—8
%迭代
jx40=jx41;
else
break;
end
end
disp(jx41');
运营结果:(n等于4,5,6时)x的转置
当n=4时:
[0.770.460.460.770.460.230.230.460.460.230.
230.460.770.460.460.77]
当n=5时:
E0.790.500.580.500.790.500.370.000.370.5
00.580.000.160.000.580.500.370.000.370.500.79
0.500.580.500.79]
当n=6时的结果:
[0.840.430.830.830.430.840.43
0.670.260.260.670.43
0.830.260.100.100.260.830.830.26
0.100.100.260.83
0.430.670.260.260.670.430.840.4
30.830.830.430.84]
(3)Gauss-Seide1迭代法代码(只列出n等于4时的代码,n等于5和6
时类似)
%取初始向量为0
gsx40=zeros(16,1);
gsx41=zeros(16,1);
gs—dx4=zeros(16,1);
while1
fori=1:16
sum=0;
forj=l:i-1
sum=sum+A4(i,j)*gsx40(j);
end
forj=i:16
sum=sum+A4(i,j)*gsx40(j);
end
gs_dx4(i)=(b4(i)-sum)/A4(i,i);
gsx41(i)=gsx40(i)+gs_dx4(i);
end
ifnorm(gsx41-gsx40)>le—8
%迭代
gsx40=gsx41;
else
break;
end
end
disp(gsx41');
运营结果:(n等于4,5,6时)x的转置
当n=4时的结果:
[0.770.460.460.77
0.460.230.230.46
0.460.230.230.46
0.770.460.460.77]
当n=5时的结果:
[0.790.500.
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 建筑预应力管桩劳务分包协议
- 湖北省襄阳市保康县第一实验小学2024-2025学年六年级上学期12月月考数学试题(含答案)
- 房屋买卖定金合同的支付地点
- 英语必修一外研版课本课件精讲攻略
- 外研版初中英语七年级上册教案总述
- 招标文件审查记录表
- 医疗机构医废管理培训
- 国防生训练协议模板
- 2023年特殊工种防突工考试题库(含答案)
- 九年级上外研版数学竞赛辅导
- 2014新农人报告-阿里研究院
- 四措两案规范标准材料模板
- 人教版新版高一英语-第一册-welcome-unitPPT课件
- 小儿中医体质辨识量表
- 各种建筑物的冷热负荷指标
- 日照钢铁 6-150028.65吨 质量证明书
- 数字经济背景下提升中华优秀传统文化影响力的思考
- 安全工程—英语双专业(双学位)培养计划(精)
- 财神正朝科仪
- 体格检查基本规范
- 毕业论文打印机皮带驱动系统能控能观和稳定性分析
评论
0/150
提交评论