数值分析实验报告2.docx
- 文档编号:6997583
- 上传时间:2023-01-15
- 格式:DOCX
- 页数:16
- 大小:33.07KB
数值分析实验报告2.docx
《数值分析实验报告2.docx》由会员分享,可在线阅读,更多相关《数值分析实验报告2.docx(16页珍藏版)》请在冰豆网上搜索。
数值分析实验报告2
数值分析上机实验报告
(2)
(注:
本实验报告中所有程序均为MATLAB语言程序)
班级:
姓名:
学号:
六章
2、用比例求根法求
在区间(1,
/2)内的一个根,直到近似根
满足精度|
(
)|<0.00001时终止运算。
程序:
clear
formatlong
x=1;
f=inline('1-x*sin(x)');
a=1;
b=pi/2;
whileabs(a-b)>0.00001
iff((a+b)/2)<0
b=(a+b)/2;
x=b;
end
iff((a+b)/2)>0
a=(a+b)/2;
x=a;
end
iff((a+b)/2)==0
x=(a+b)/2;
break
end
end
x
运行结果:
x=
1.114153168596455
4、比较一下两种求
+10
-2=0的根到三位有效数字所需的计算量。
(1)在区间(0,1)内用二分法;
(2)用迭代法
=(2-
)/10,取初值
=0。
(1)程序:
formatlong
a=0;
b=1;
R=1;
k=0;
f=inline('exp(x)+10*x-2');
whileR>1e-5
c=(a+b)/2;
iff(a)*f(c)>0;
a=c;
else
b=c;
end
R=b-a;
k=k+1;
end
k
x=c
运行结果:
k=
17
x=
0.09052276611328
首先初步估计精确值在0.09左右,故取误差限选择为10^-5,,得到的计算次数为17次。
(2)程序:
formatlong
f=inline('(2-exp(x))/10');
disp('x=');
x=feval(f,0);
disp(x);
Eps=1E-5;
i=1;
while1
x0=x;
i=i+1;
x=feval(f,x);
disp(x);
ifx>1E10;
break;
end;
ifabs(x-x0) break; end end i,x 输出结果: x= 0.10000000000000 0.08948290819244 .0906******** 0.09051261667437 0.09052646805264 0.09052495168284 i= 6 x= 0.09052495168284 若要保留三位有效数字,则误差限应选择为10^-5,得到的迭代次数为5次。 7、用下列方法求 =0在 =2附近的根,根的准确值 =1.87938524…,要求计算结果准确到四位有效数字。 (1)用newton法; (2)用弦截法,取 =2, =1.9; (3)用抛物线法,取 =1, =3, =2。 (1)Newton法程序: clear m=input('m='); x (1)=2; f=inline('x^3-3*x-1'); f1=inline('3*x^2-3'); forn=1: m x(n+1)=x(n)-f(x(n))/f1(x(n)); end x(n+1) 运行结果: m=10 ans= 1.879385241571817 (2)弦截法程序: clear m=input('m='); x (1)=2; x (2)=1.9; f=inline('x^3-3*x-1'); forn=2: m x(n+1)=x(n)-(f(x(n))/(f(x(n))-f(x(n-1))))*(x(n)-x(n-1)); end x(n+1) 运行结果: m=9 ans= 1.879385241571817 (3)抛物线法程序: 首先建立M文件: functionw=W(a,b) f=inline('x^3-3*x-1'); w=(f(b)-f(a))/(b-a); 再有程序: clear m=input('m='); x (1)=1; x (2)=3; x(3)=2; f=inline('x^3-3*x-1'); forn=3: m q(n)=W(x(n),x(n-1)); p(n)=(W(x(n),x(n-2))-W(x(n),x(n-1)))/(x(n-2)-x(n-1)); w(n)=q(n)+p(n)*(x(n)-x(n-1)); s(n)=sqrt(w(n)^2-4*f(x(n))*p(n)); ifw(n)+s(n)>0 x(n+1)=x(n)-2*f(x(n))/(w(n)+sqrt(w(n)^2-4*f(x(n))*p(n))); end ifw(n)-s(n)>0 x(n+1)=x(n)-2*f(x(n))/(w(n)-sqrt(w(n)^2-4*f(x(n))*p(n))); end end x=x(n+1) 运行结果: m=10 x= 1.879385241570596 13.应用Newton法于方程 导出求 的迭代公式,并求出 的值 程序: clear m=input('m='); x (1)=11; f=inline('1-115/x^2'); f1=inline('2*115*x^(-3)'); forn=1: m x(n+1)=x(n)-f(x(n))/f1(x(n)); end sqrta=x(n+1) 运行结果: m=10 sqrta= 10.723805294763608 七、八、九章 求解下面线性方程的解: A是一个主对角线为4,主对角线上下各三条对角线为1的100阶的方阵,b是一个1到100的一列矩阵。 利用不同种方法就解线性方程组。 Lu分解程序: clear formatshort n=input('n='); a=ones(n); b=diag(diag(a,-1),-1); c=diag(diag(a,-2),-2); d=diag(diag(a,1),1); e=diag(diag(a,2),2); g=diag(diag(a,3),3); h=diag(diag(a,-3),-3); f=4*eye(n); A=b+c+d+e+f+g+h; b=[1: 100]'; [L,U]=lu(A); x=U\(L\b); x' 运行结果: n=100 ans= Columns1through11 0.01700.18550.32260.42390.49450.58860.69640.80500.90330.99931.0978 Columns12through22 1.19961.30081.40061.49981.59971.69991.80011.90012.00002.09992.2000 Columns23through33 2.30002.40002.50002.60002.70002.80002.90003.00003.10003.20003.3000 Columns34through44 3.40003.50003.60003.70003.80003.90004.00004.10004.20004.30004.4000 Columns45through55 4.50004.60004.70004.80004.90005.00005.10005.20005.30005.40005.5000 Columns56through66 5.60005.70005.80005.90006.00006.10006.20006.30006.40006.50006.6000 Columns67through77 6.69996.80006.90017.00027.09997.19967.29987.40067.50097.59957.6976 Columns78through88 7.79907.90398.00548.09648.18638.29498.42488.52928.57808.61588.7831 Columns89through99 9.04009.18158.92308.75859.153710.391110.27678.85456.469110.266214.0251 Column100 17.3099 Cholesky分解: clear formatshort n=input('n='); a=ones(n); b=diag(diag(a,-1),-1); c=diag(diag(a,-2),-2); d=diag(diag(a,1),1); e=diag(diag(a,2),2); g=diag(diag(a,3),3); h=diag(diag(a,-3),-3); f=4*eye(n); A=b+c+d+e+f+g+h; b=[1: 100]'; R=chol(A); x=R\(R'\b); x' 运行结果: n=100 ans= Columns1through11 0.01700.18550.32260.42390.49450.58860.69640.80500.90330.99931.0978 Columns12through22 1.19961.30081.40061.49981.59971.69991.80011.90012.00002.09992.2000 Columns23through33 2.30002.40002.50002.60002.70002.80002.90003.00003.10003.20003.3000 Columns34through44 3.40003.50003.60003.70003.80003.90004.00004.10004.20004.30004.4000 Columns45through55 4.50004.60004.70004.80004.90005.00005.10005.20005.30005.40005.5000 Columns56through66 5.60005.70005.80005.90006.00006.10006.20006.30006.40006.50006.6000 Columns67through77 6.69996.80006.90017.00027.09997.19967.29987.40067.50097.59957.6976 Columns78through88 7.79907.90398.00548.09648.18638.29498.42488.52928.57808.61588.7831 Columns89through99 9.04009.18158.92308.75859.153710.391110.27678.85456.469110.266214.0251 Column100 17.3099 Jacobi方法迭代 函数M文件: function[y,n]=jacobi(A,b,x0,eps) ifnargin==3 eps=1.0e-6; elseifnargin<3 error return end D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=D\(L+U); f=D\b; y=B*x0+f; n=1; whilenorm(y-x0)>=eps x0=y; y=B*x0+f; n=n+1; end 程序: clear formatshort n=input('n='); a=ones(n); b=diag(diag(a,-1),-1); c=diag(diag(a,-2),-2); d=diag(diag(a,1),1); e=diag(diag(a,2),2); g=diag(diag(a,3),3); h=diag(diag(a,-3),-3); f=4*eye(n); A=b+c+d+e+f+g+h; b=[1: n]'; m=diag(a); [x,n]=jacobi(A,b,m,1.0e-3); x' n 运行结果: 因该方程的系数矩阵A不满足相关条件,故Jacobi方法迭代对该方程不适用 Gauss-Serdel迭代方法 函数M文件: function[y,n]=gauseidel(A,b,x0,eps) ifnargin==3 eps=1.0e-6; elseifnargin<3 error return end D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); G=(D-L)\U; f=(D-L)\b; y=G*x0+f; n=1; whilenorm(y-x0)>=eps x0=y; y=G*x0+f; n=n+1; end 命令程序: clear formatlong n=input('n='); a=ones(n); b=diag(diag(a,-1),-1); c=diag(diag(a,-2),-2); d=diag(diag(a,1),1); e=diag(diag(a,2),2); g=diag(diag(a,3),3); h=diag(diag(a,-3),-3); f=4*eye(n); A=b+c+d+e+f+g+h; b=[1: n]'; m=diag(a); [x,n]=gauseidel(A,b,m,1.0e-10); x' n 运行结果: n=100 ans= Columns1through10 0.01700.18550.32260.42390.49450.58860.69640.80500.90330.9993 Columns11through20 1.09781.19961.30081.40061.49981.59971.69991.80011.90012.0000 Columns21through30 2.09992.20002.30002.40002.50002.60002.70002.80002.90003.0000 Columns31through40 3.10003.20003.30003.40003.50003.60003.70003.80003.90004.0000 Columns41through50 4.10004.20004.30004.40004.50004.60004.70004.80004.90005.0000 Columns51through60 5.10005.20005.30005.40005.50005.60005.70005.80005.90006.0000 Columns61through70 6.10006.20006.30006.40006.50006.60006.69996.80006.90017.0002 Columns71through80 7.09997.19967.29987.40067.50097.59957.69767.79907.90398.0054 Columns81through90 8.09648.18638.29498.42488.52928.57808.61588.78319.04009.1815 Columns91through100 8.92308.75859.153710.391110.27678.85456.469110.266214.025117.3099 n= 46 幂法: 建立matlab的函数文件eig_power.m如下: function[V,D]=eig_power(A) Maxtime=100; Eps=1E-5; n=length(A); V=ones(n,1); k=0; m0=0; whilek<=Maxtime v=A*V; [vmax,i]=max(abs(v)); m=v(i); V=v/m; ifabs(m-m0) break; end m0=m; k=k+1; end D=m; 在命令窗口中调用函数文件eig_power.m: a=linspace(4,4,100); c=linspace(1,1,99); d=linspace(1,1,98); e=linspace(1,1,97); f=diag(a); g=diag(c,1); h=diag(c,-1); i=diag(d,2); j=diag(d,-2); k=diag(e,3); l=diag(e,-3); A=f+g+h+i+j+k+l; [V,D]=eig_power(A); V' D 输出结果: ans= Columns1through9 0.55000.68000.81000.94000.97000.99001.00001.00001.0000 Columns10through18 1.00001.00001.00001.00001.00001.00001.00001.00001.0000 Columns19through27 1.00001.00001.00001.00001.00001.00001.00001.00001.0000 Columns28through36 1.00001.00001.00001.00001.00001.00001.00001.00001.0000 Columns37through45 1.00001.00001.00001.00001.00001.00001.00001.00001.0000 Columns46through54 1.00001.00001.00001.00001.00001.00001.00001.00001.0000 Columns55through63 1.00001.00001.00001.00001.00001.00001.00001.00001.0000 Columns64through72 1.00001.00001.00001.00001.00001.00001.00001.00001.0000 Columns73through81 1.00001.00001.00001.00001.00001.00001.00001.00001.0000 Columns82through90 1.00001.00001.00001.00001.00001.00001.00001.00001.0000 Columns91through99 1.00001.00001.00001.00000.99000.97000.94000.81000.6800 Column100 0.5500 D= 10 其中D是最大特征值,V则是特征向量。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 实验 报告
![提示](https://static.bdocx.com/images/bang_tan.gif)