东北大学秦皇岛分校数值分析实验报告.docx
- 文档编号:24685878
- 上传时间:2023-05-31
- 格式:DOCX
- 页数:23
- 大小:716.68KB
东北大学秦皇岛分校数值分析实验报告.docx
《东北大学秦皇岛分校数值分析实验报告.docx》由会员分享,可在线阅读,更多相关《东北大学秦皇岛分校数值分析实验报告.docx(23页珍藏版)》请在冰豆网上搜索。
东北大学秦皇岛分校数值分析实验报告
数值分析
实验报告
学号:
[2133326]
姓名:
[孙丽]
提交日期:
[2015-6-25]
成绩:
东北大学秦皇岛分校
【实验编号】
1
【实验名称】
多项式插值
【实验目的】
1)掌握拉格朗日插值和牛顿插值;
2)掌握用mathematica编写多项式插值的应用程序;
【实验内容】
给定函数
,将区间[-3,3]剖分为12等份,取节点
,分别构造拉格朗日插值和牛顿插值多项式
和
,比较两个多项式的计算时间和精度,最后绘出
和
的图形。
【实验结果】
1.f[x_]:
=Sqrt[2Pi]Exp[-x^2/2]
g=Plot[f[x],{x,-3,3},PlotStyle->{Red,Thick},BaseStyle->FontSize->15,
AxesLabel->{x,y},AxesStyle->Thick]
2.lagrange[x_,xp_,yp_]:
=Module[{n,s=0},n=Length[xp];
Do[s=s+Product[(x-xp[[j]])/(xp[[i]]-xp[[j]]),{j,i-1}]
Product[(x-xp[[j]])/(xp[[i]]-xp[[j]]),{j,i+1,n}]yp[[i]],{i,n}];s]
newton[x_,xp_,yp_]:
=Module[{n,p,yp0=yp},n=Length[xp];
Do[Do[yp0[[j]]=(yp0[[j]]-yp0[[j-1]])/(xp[[j]]-xp[[j-i]]),
{j,n,i+1,-1}],{i,n-1}];p=yp0[[n]];Do[p=p(x-xp[[i]])+yp0[[i]],{i,n-1,1,-1}];p]
xd=Range[-3.,3.,0.5];yd=Map[f,xd];dat=Transpose[{xd,yd}];
n=Length[dat]
l[x_]=lagrange[x,xd,yd];
Clear[p]
p[x_]=newton[x,xd,yd];
Plot[l[x],{x,xd[[1]],xd[[n]]}]//Timing
3.Plot[p[x],{x,xd[[1]],xd[[n]]}]//Timing
4.h[x_]=InterpolatingPolynomial[dat,x];
Plot[h[x],{x,xd[[1]],xd[[n]]}]//Timing
5.g=Plot[l[x]-p[x],{x,xd[[1]],xd[[n]]},PlotRangeAll,
PlotStyleBlue,BaseStyle{Bold,FontSize14},
AxesStyleThick,AxesLabel{x,error}]
6.Plot[h[x]-p[x],{x,xd[[1]],xd[[n]]},PlotRangeAll]
7.Clear[p1]
p1[x_]=p[x]
8.Clear[p]
g=Plot[p1[x],{x,xd[[1]],xd[[n]]},PlotStyle{Red,Thick},
BaseStyleFontSize15,AxesStyleThick,
AxesLabel{x,p[x]},
Epilog{PointSize[Large],Blue,Point/2@dat}]
9.error=.
g=Plot[f[x]-p1[x],{x,xd[[1]],xd[[n]]},PlotRangeAll,
PlotStyle{Red,Thick},BaseStyle{Bold,FontSize16},
AxesLabel{x,error},AxesStyleThick]
【实验编号】
2
【实验名称】
数值积分
【实验目的】
1)掌握复合梯形公式和复合辛普森公式;
2)掌握用mathematica编写数值积分的应用程序;
【实验内容】
给定函数
,动画演示复合梯形公式和符合辛普森公式。
【实验结果】
1.ntos[n_Integer]:
=Module[{s},s=IntegerDigits[n];
If[n<10,"0"<>FromCharacterCode[s+48],FromCharacterCode[s+48]]]
ntos[45]
2.f[x_]:
=12+4Sin[2x]Exp[x/10]
t1=Plot[f[x],{x,0,11},AxesLabel->{x,y},
PlotRange->{{0,12},{0,25}},
PlotStyle->{Thickness[0.01],RGBColor[1,0,1]},AxesStyle->Thick]
3.a=0;b=11;
Manipulate[dat=Table[{xp=a+((b-a)i)/n,f[xp]},{i,0,n}];
l1=Interpolation[dat,InterpolationOrder->1];
g1=Plot[l1[x],{x,a,b},PlotRange->{{0,12},{0,25}},
Filling->Axis,FillingStyle->Green];
Show[t1,g1,PlotLabelntos[n+1]<>"PointsFigure"],{n,2,40,2}]
4.ntos[n_Integer]:
=Module{s},s=IntegerDigits[n];
If[n<10,"0"<>FromCharacterCode[s+48],FromCharacterCode[s+48]]]
g[x_]:
=12+4Sin[2x]Exp[x/10]
t2=Plot[g[x],{x,0,11},AxesLabel{x,g[x]},
PlotRange{{0,12},{0,25}},
PlotStyle{Thickness[0.01],RGBColor[1,0,1]},AxesStyleThick,
BaseStyleFontSize16]
a=0;
b=11;
Manipulate[dat=Table[{xp=a+((b-a)i)/n,f[xp]},{i,0,n}];
l2=Interpolation[dat,InterpolationOrder2];
g2=Plot[l2[x],{x,a,b},PlotRange{{0,12},{0,25}},
FillingAxis,FillingStyleGreen];
Show[t2,g2,PlotLabelntos[n+1]<>"PointsFigure"],{n,2,30,2}]
【实验编号】
3
【实验名称】
线性方程组求解
【实验目的】
1)掌握高斯类消去法求解线性方程组;
2)掌握用mathematica编写直接求解线性方程组的应用程序;
【实验内容】
列主元高斯消去法求解线性代数方程组Ax=b,其中
,
。
要求:
(1)对矩阵A,找出第2列中从第2个元素开始的绝对值最大元素所在位置,并与第2行进行交换。
(2)使用列主元高斯消去法将矩阵A化为上三角形,矩阵A下三角部分存放消元过程中的变换元素。
(3)回代求出线性方程组的解,存放在向量b中。
(4)将求解Ax=b的列主元高斯消去法写成函数形式。
【实验结果】
(1)A={{2,1,-3,-1},{3,1,0,7},{-1,2,4,-2},{1,0,-1,5}}
i=2;m=First[Dimensions[A]];
temp=Flatten[Abs[Take[A,{i,m},{i}]]];
imax=Position[temp,Max[temp]][[1,1]]+i-1;
{A[[i]],A[[imax]]}={A[[imax]],A[[i]]}
(2)A={{2,1,-3,-1},{3,1,0,7},{-1,2,4,-2},{1,0,-1,5}};
b={1,2,-1,5};
{m,n}=Dimensions[A];
Do[temp=Flatten[Abs[Take[A,{k,m},{k}]]];
kmax=Position[temp,Max[temp]][[1,1]]+k-1;
If[k!
=kmax,{A[[k]],A[[kmax]]}={A[[kmax]],A[[k]]};
{b[[k]],b[[kmax]]}={b[[kmax]],b[[k]]}];
If[A[[k,k]]!
=0,Do[A[[i,k]]=A[[i,k]]/A[[k,k]];
Do[A[[i,j]]=A[[i,j]]-A[[i,k]]*A[[k,j]],{j,k+1,n}];
b[[i]]=b[[i]]-A[[i,k]]*b[[k]],{i,k+1,m}],
Print["Error"];Break[]],{k,m-1}];
(3)b[[m]]=b[[m]]/A[[m,m]];
Do[b[[i]]=(b[[i]]-Sum[A[[i,j]]*b[[j]],{j,i+1,m}])/A[[i,i]],{i,m-1,1,-1}];b
(4)gausspiv[A_,b_]:
=Module[{m,n,temp,kmax,A0=A,b0=b},
{m,n}=Dimensions[A0];
Do[temp=Flatten[Abs[Take[A0,{k,m},{k}]]];
kmax=Position[temp,Max[temp]][[1,1]]+k-1;
If[kkmax,{A0[[k]],A0[[kmax]]}={A0[[kmax]],A0[[k]]};
{b0[[k]],b0[[kmax]]}={b0[[kmax]],b0[[k]]}];
If[A0[[k,k]]0,Do[A0[[i,k]]=A0[[i,k]]/A0[[k,k]];
Do[A0[[i,j]]=A0[[i,j]]-A0[[i,k]]*A0[[k,j]],{j,k+1,n}];
b0[[i]]=b0[[i]]-A0[[i,k]]*b0[[k]],{i,k+1,m}],
Print["Error"];Break[]],{k,m-1}];
b0[[m]]=b0[[m]]/A0[[m,m]];
Do[b0[[i]]=(b0[[i]]-Sum[A0[[i,j]]*b0[[j]],{j,i+1,m}])/A0[[i,i]],{i,m-1,1,-1}];
b0]
x=gausspiv[A,b]
【实验编号】
4
【实验名称】
非线性方程求根
【实验目的】
1)掌握二分法和牛顿法求解非线性方程的根;
2)掌握用mathematica编写非线性方程求根的应用程序;
【实验内容】
给定非线性方程
,完成以下工作:
(1)绘出函数
的图形,确定有根区间。
(2)使用二分法求出方程的根,要求精度
(3)使用牛顿法求出方程的根,要求精度
【实验结果】
(1)bisection[f_,a_,b_,eps_,kmax_]:
=Module[{k=0,v=1,a0=a,b0=b,va,vb,c},While[Sign[f[a]]Sign[f[b]],a=Input["Inputa:
"];
b=Input["Inputb:
"]];
Print["a=",a,"b=",b];
va=f[a0];
vb=f[b0];
While[Abs[v]eps&&k v=f[c]; If[Sign[va]Sign[v],a0=c; va=v,b0=c; vb=v]; k++; Print["k=",k,"c=",c,"f(",c,")=",v]]; {c,v}] f[x_]: =Sin[10x]-xExp[-x]-x; res=bisection[f,0.1,0.5,10^-8,100] (2) newton[f_,x0_,eps_,kmax_]: =Module[{k=0,v=1,x00=x0,x1,fd,t},fd[t_]=D[f[t],t];v=f[x00]//N; Print["k=",k,"x0=",x00,"f(",x00,")=",v]; While[Abs[v]eps&&k v=f[x1]; x00=x1; k++; Print["k=",k,"x1=",x1,"f(",x1,")=",v]]; If[Abs[v]>eps,Print["Newton'smethoddidnotproduceanaccurateapproximation"]]; {x1,v}] f[x_]: =Sin[10x]-xExp[-x]-x; res=newton[f,0.4,10^-8,100] (3) res=newton[f,-1,10^-8,100] (4) res=newton[f,-0.1,10^-8,100] (5)res=newton[f,-0.5,10^-8,20]
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 东北大学 秦皇岛 分校 数值 分析 实验 报告
![提示](https://static.bdocx.com/images/bang_tan.gif)