计算方法实验.docx
- 文档编号:9140824
- 上传时间:2023-02-03
- 格式:DOCX
- 页数:17
- 大小:31.57KB
计算方法实验.docx
《计算方法实验.docx》由会员分享,可在线阅读,更多相关《计算方法实验.docx(17页珍藏版)》请在冰豆网上搜索。
计算方法实验
计算方法实验报告
目录
计算方法实验报告-1-
实验一舍入误差与数值稳定性-2-
实验二方程求根-3-
实验三线性方程组数值解法-4-
实验四插值法-8-
实验五曲线拟合-9-
实验六数值积分与数值微分-11-
实验一舍入误差与数值稳定性
设
=
已知其精确值为
(
.
1)编制从大到小的顺序计算
的程序;
2)编制按从小到大的顺序计算
的程序;
3)按2种顺序分别计算
,
,
,并指出有效位数。
思路分析:
两种不同的累加方式。
解:
1)程序如下:
(预设N=100)
N=100;
S=0;
while(N-1)
S=S+1/(N^2-1);
N=N-1;
end
fprintf('%f',S);
2)程序如下:
N=100;
S=0;
forj=2:
N
S=S+1/(j^2-1);
end
fprintf('%f',S);
3)降序:
=0.749000
=0.749900
=0.749967
升序:
=0.749000
=0.749900
=0.749967
实验二方程求根
用牛顿法求下列方程的根:
思路分析:
首先应是牛顿迭代式,之后判断收敛即可。
解:
程序如下:
eps=5e-6;
delta=1e-6;
N=100;
k=0;
x0=1;
while
(1)
x1=x0-func1(x0)/func2(x0);
k=k+1;
if(k>N|abs(x1) disp('newtonmethodfailed') break; end ifabs(x1)<1 d=x1-x0; else d=(x1-x0)/x1; end x0=x1; if(abs(d) break; end end fprintf('%f',x0); functiony=func2(x) y=2*x-exp(x); end functiony=func1(x) y=x*x-exp(x); end 运行结果: -0.703467 实验三线性方程组数值解法 1.编写用追赶法解三对角线性方程组的程序,并解下列方程: Ax=b,其中, =[-41 =[-27 1-41-15 1-41-15 ………… 1-41-15 1-4]-15] 思路分析: 按照追赶法的递推公式求解。 解: 程序如下: functionx=Trid(a,b,c,d) %追赶法求解三对角的线性方程组Ax=d %b为主对角线元素,a,c分别为次对角线元素,d为右端项 %A=[b1c1 %a2b2c2 %...... %a_(n-1)b_(n-1)c_(n-1) %a_(n)b_(n)] %b=[b1...b_(n)] %a=[0a2...a_(n)] %c=[c1...c_(n-1)] n=length(b); u (1)=b (1); fori=2: n l(i)=a(i)/u(i-1); u(i)=b(i)-l(i)*c(i-1); end y (1)=d (1); fori=2: n y(i)=d(i)-l(i)*y(i-1); end x(n)=y(n)/u(n); fori=n-1: -1: 1 x(i)=(y(i)-c(i)*x(i+1))/u(i); end fori=1: n fprintf('x[%d]=%f\n',i,x(i)); end 运行结果: x[1]=8.705758,x[2]=7.823032,x[3]=7.586371,x[4]=7.522453,x[5]=7.503440,x[6]=7.491306,x[7]=7.461785,x[8]=7.355835,x[9]=6.961556,x[10]=5.490389 2.分别用雅克比迭代法与高斯-赛德尔迭代法解下列方程组: RI=V,其中 R=[31-13000-10000V=[-15 -1335-90-11000027 0-931-1000000-23 00-1079-30000-90 000-3057-70-50-20 0000747-300012 00000-304100-7 0000-50027-27 0000000-229]-10] 思路分析: 分别按照雅克比与高斯迭代法递推式编制程序。 解: 程序如下: 雅克比迭代法 functionjacobi(A,b,max,eps) n=length(A);x=zeros(n,1);x1=zeros(n,1);k=0; while1 x1 (1)=(b (1)-A(1,2: n)*x(2: n,1))/A(1,1); fori=2: n-1 x1(i)=(b(i)-A(i,1: i-1)*x(1: i-1,1)-A(i,i+1: n)*x(i+1: n,1))/A(i,i); end x1(n)=(b(n)-A(n,1: n-1)*x(1: n-1,1))/A(n,n); k=k+1; ifsum(abs(x1-x)) fprint('number=%d\n',k); break; end ifk>=max fprintf('themethodisdisconvergent\n'); break; end x=x1; end ifk fori=1: n fprintf('x[%d]=%f\n',i,x1(i)); end end 运行结果: x[1]=-0.200550,x[2]=0.368393,x[3]=-0.731860,x[4]=-0.300318,x[5]=-0.446577,x[6]=0.399384,x[7]=0.121500,x[8]=0.151792,x[9]=-0.334359 高斯-赛德尔迭代法 functionGauseSeidel(A,b,max,eps) n=length(A);x=zeros(n,1);x1=zeros(n,1);k=0; while1 x1 (1)=(b (1)-A(1,2: n)*x(2: n,1))/A(1,1); fori=2: n-1 x1(i)=(b(i)-A(i,1: i-1)*x1(1: i-1,1)-A(i,i+1: n)*x(i+1: n,1))/A(I,I); end x1(n)=(b(n)-A(n,1: n-1)*x1(1: n-1,1))/A(n,n); k=k+1; ifsum(abs(x1-x)) fprintf('number=%d\n',k); break; end ifk>=max fprintf('themethodisdiscomvergent\n'); break; end x=x1; end ifk fori=1: n fprint('x[%d]=%f\n',i,x1(i)); end end 运行结果: x[1]=-0.200550,x[2]=0.368393,x[3]=-0.731860,x[4]=-0.300318,x[5]=-0.446577,x[6]=0.399384,x[7]=0.121500,x[8]=0.151792,x[9]=-0.334359 实验四插值法 按下列数据做五次插值,并求x1=0.46,x2=0.55,x3=0.60时的函数近似值。 Xi 0.30 0.42 0.50 0.58 0.66 0.72 yi 1.04403 1.08462 1.11803 1.15603 1.19817 1.23223 思路分析: 用牛顿插值法多项式来解决此问题。 解: 程序如下: functionnewtint(x,y,xhat) n=length(y); c=y(: ); forj=2: n fori=n: -1: j c(i)=(c(i)-c(i-1))/(x(i)-x(i-j+1)); end end yhat=c(n); fori=n-1: -1: 1 yhat=yhat*(xhat-x(i))+c(i); end fprintf('N(%f)=%f',xhat,yhat); 运行结果: x=0.46,y=1.100742; x=0.55,y=1.141271; x=0.60,y=1.166194; 实验五曲线拟合 试分别用抛物线y=a+bx+c 和指数曲线y=a 拟合下列数据: X 1 1.5 2 2.5 3 3.5 4 4.5 y 33.4 79.5 122.65 159.05 189.15 214.15 238.65 252.50 x 5 5.5 6 6.5 7 7.5 8 y 267.55 280.5 296.65 301.40 310.4 318.15 325.15 比较2个拟合函数的优劣。 解: 1.抛物线拟合程序如下: functionZXE(x,y,m) S=zeros(1,2*m+1);T=zeros(m+1,1); fork=1: 2*m+1 S(k)=sum(x.^(k-1)); end fork=1: m+1 T(k)=sum(x.^(k-1).*y); end A=zeros(m+1,m+1);a=zeros(m+1,1); fori=1: m+1 forj=1: m+1 A(i,j)=S(i+j-1); end end a=A\T; fork=1: m+1 fprintf('a[%d]=%f\n',k,a(k)); end 运行结果: a[1]=-45.333297,a[2]=94.230200,a[3]=-6.131610;即: y=-45.333297+94.230200x-6.131610 2.指数曲线拟合程序如下: functionZXEI(x,y,m) y=log(y); S=zeros(1,2*m-1);T=zeros(m,1); fork=1: 2*m-1 S(k)=sum(x.^(k-1)); end fork=1: m T(k)=sum(x.^(k-1).*y); end A=zeros(m,m);a=zeros(m,1); fori=1: m forj=1: m A(i,j)=S(i+j-1); end end a=A\T; a (1)=exp(a (1)); fork=1: m fprintf('a[%d]=%f\n',k,a(k)); end 运行结果: a[1]=67.402,a[2]=0.238960,即: y=67.402 实验六数值积分与数值微分 用复化梯形公式和复化辛普森公式计算积分 和 观察n为多少时所得近似值具有6位有效数字。 解: MATLAB程序如下: 复化梯形公式 functiontrap(a,b) n=1030; fori=1: 3 x=a: (b-a)/n: b; h=(b-a)/n; s=(h/2)*(f(a)+2*sum(f(x(2: 1: n-1)))+f(b)); fprintf('s(%d)=%f\n',n,s); n=n+1; end functiony=f(x) y=(1+cos(x).^2).^0.5; 运行结果: s(1030)=1.908574 s(1031)=1.908575 s(1032)=1.908577 n=1030时才有6位有效数字。 复化辛普森公式 functionsimpson(a,b) n=2; fori=1: 3 x=a: (b-a)/(2*n): b; m=2*n+1; h=(b-a)/n; s=(h/6)*(f(a)+2*sum(f(x(3: 2: m-2)))+4*sum(f(x(2: 2: m-1)))+f(b)); fprintf('s(%d)=%f\n',n,s); n=n*2; end functiony=f(x) y=(1+cos(x).^2).^0.5; 运行结果: s (2)=1.910141 s(4)=1.910099 s(8)=1.910099 n=4时即有6位有效数字; 复化梯形公式 functiontrap(a,b) n=1954; fori=1: 3 x=a: (b-a)/n: b; h=(b-a)/n; s=(h/2)*(1+2*sum(f(x(2: 1: n-1)))+f(b)); fprintf('s(%d)=%f\n',n,s); n=n+1; end functiony=f(x) y=tan(x)/x; 运行结果: s(1954)=0.848456 s(1955)=0.848456 s(1956)=0.848456 n=1954时才有6位有效数字。 复化辛普森公式 functionsimpson(a,b) n=4; fori=1: 3 x=a: (b-a)/(2*n): b; m=2*n+1; h=(b-a)/n; s=(h/6)*(1+2*sum(f(x(3: 2: m-2)))+4*sum(f(x(2: 2: m-1)))+f(b)); fprintf('s(%d)=%f\n',n,s); n=n*2; end functiony=f(x) y=tan(x)./x; 运行结果: s(4)=0.848972 s(8)=0.848967 s(16)=0.848967 n=8时即有6位有效数字;
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 实验