数值方法计算实习题信计101班和数应1.docx
- 文档编号:11981459
- 上传时间:2023-04-16
- 格式:DOCX
- 页数:26
- 大小:87.17KB
数值方法计算实习题信计101班和数应1.docx
《数值方法计算实习题信计101班和数应1.docx》由会员分享,可在线阅读,更多相关《数值方法计算实习题信计101班和数应1.docx(26页珍藏版)》请在冰豆网上搜索。
数值方法计算实习题信计101班和数应1
信计091
要求:
1、用Matlab语言或你熟悉的其他算法语言编写程序,使之尽可能具有通用性;2、根据上机计算实践,对所使用的数值方法的特点、性质、有效性、误差和收敛性等方面进行必要的讨论和分析;3、完成计算后写出实验报告,内容包括:
课题名称、解决的问题、采用的数值方法、算法程序、数值结果、对实验结果的讨论和分析等;4、特别说明:
严禁抄袭,否则一经发现,所有雷同实验报告最多评为及格。
一、下表给出了飞行中鸭子的上部形状的节点数据,试用三次样条插值函数(自然边界条件)和20次Lagrange插值多项式对数据进行插值。
用图示出给定的数据,以及
和
。
0.9
1.3
1.9
2.1
2.6
3.0
3.9
4.4
4.7
5.0
6.0
1.3
1.5
1.85
2.1
2.6
2.7
2.4
2.15
2.05
2.1
2.25
7.0
8.0
9.2
10.5
11.3
11.6
12
12.6
13.0
13.3
2.3
2.25
1.95
1.4
0.9
0.7
0.6
0.5
0.4
0.25
解:
>>x=[0.91.31.92.12.63.03.94.44.75.06.07.08.09.210.511.311.61212.613.013.3];
>>y=[1.31.51.852.12.62.72.42.152.052.12.252.32.251.951.40.90.70.60.50.40.25];
(1)三次样条插值法
在MATLAB中编写m文件
function[f,f0]=scyt(x,y,y2_1,y2_N,x0)
%y2_1和y2_N分别为自然边界条件
%插值点x的坐标:
x0
symst;
f=0.0;f0=0.0;
if(length(x)==length(y))
n=length(x);
else
disp('x和y的维数不相等');
return;
end
fori=1:
n
if(x(i)<=x0)&&(x(i+1)>=x0)
index=i;
break;
end
end
A=diag(2*ones(1,n));
A(1,2)=1;
A(n,n-1)=1;
u=zeros(n-2,1);
lamda=zeros(n-1,1);
c=zeros(n,1);
fori=2:
n-1
u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));
lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));
c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));
A(i,i+1)=u(i-1);
A(i,i-1)=lamda(i);
end
c
(1)=3*(y
(2)-y
(1))/(x
(2)-x
(1))-(x
(2)-x
(1))*y2_1/2;
c(n)=3*(y(n)-y(n-1))/(x(n)-x(n-1))-(x(n)-x(n-1))*y2_N/2;
m=zgf(A,c);
h=x(index+1)-x(index);
f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h;
f0=subs(f,'t',x0);
其中的zgf函数为追赶法,其程序为
functionx=zgf(A,b)
n=rank(A);
for(i=1:
n)
if(A(i,i)==0)
disp('Error:
对角有元素为0!
');
return;
end
end;
d=ones(n,1);
a=ones(n-1,1);
c=ones(n-1);
for(i=1:
n-1)
a(i,1)=A(i+1,i);
c(i,1)=A(i,i+1);
d(i,1)=A(i,i);
end
d(n,1)=A(n,n);
for(i=2:
n)
d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1))*c(i-1,1);
b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1))*b(i-1,1);
end
x(n,1)=b(n,1)/d(n,1);
for(i=(n-1):
-1:
1)
x(i,1)=(b(i,1)-c(i,1)*x(i+1,1))/d(i,1);
end
在MATLAB中输入指令
>>[f,f0]=scyt(x,y,0,0)
f=
1000/729*(27/5*t-1377/100)*(t-39/10)^2+1000/729*(522/25-24/5*t)*(t-3)^2+100/81*(-6396162352027119/288230376151711744*t+19188487056081357/288230376151711744)*(39/10-t)^2-100/81*(-176836856862157557/90071992547409920+4534278381080963/9007199254740992*t)*(t-3)^2
f0=
2.5851得三次样条插值函数
S(x)=1000/729*(27/5*x-1377/100)*(x-39/10)^2+1000/729*(522/25-24/5*x)*(x-3)^2+100/81*(-6396162352027119/288230376151711744*x+19188487056081357/288230376151711744)*(39/10-x)^2-100/81*(-176836856862157557/90071992547409920+4534278381080963/9007199254740992*x)*(x-3)^2
>>xi=0.9:
0.01:
13.3;yi=interp1(x,y,xi,'spline');
>>title('试验一--三次样条插值图示')
(2)用拉格朗日法插值
%定义Lagrange程序
functionf=Language(x,y,x0)
symst;
if(length(x)==length(y))
n=length(x);
else
disp('x和y的维数不相等');
return;
end
f=0.0;
for(i=1:
n)
l=y(i);
for(j=1:
i-1)
l=l*(t-x(j))/(x(i)-x(j));
end;
for(j=i+1:
n)
l=l*(t-x(j))/(x(i)-x(j));
end;
f=f+l;
simplify(f);
if(i==n)
if(nargin==3)
f=subs(f,'t',x0);
else
f=collect(f);
f=vpa(f,6);
end
end
end
>>Language(x,y)
ans=
52462.6*t+189995.*t^3-189851.*t^4+136778.*t^5-11.3161*t^12-.277283e-6*t^18+1.18284*t^13-73866.6*t^6+.111076e-4*t^17-.976904e-1*t^14+.427949e-8*t^19-.307453e-10*t^20+30677.6*t^7+2564.20*t^9-9968.98*t^8+.628590e-2*t^15-525.813*t^10-9652.78-.308159e-3*t^16+86.2514*t^11-128683.*t^2
二、已知Wilson矩阵
,且向量
,则方程组
有准确解
。
⑴用Matlab内部函数求
,
的所有特征值和
;
⑵令
,解方程组
,并求出向量
和
,从理论结果和实际计算结果两方面分析方程组
解的相对误差
与
的相对误差
的关系;
⑶再改变扰动矩阵
(其元素的绝对值不超过0.005),重复第2问。
解:
(1)A=[10787;7565;86109;75910];b=[32;23;33;31];
M=det(A)
M=
1
A的所以特征值:
>>D=eig(A)
D=
0.0102
0.8431
3.8581
30.2887
条件数>>n=30.2887/0.0102
n=
2.9695e+003
所以A的行列式为1,cond(A)2=2.9695e+003
(2)>>B=[1078.17.2;7.085.0465;85.989.899;6.994.9999.98];
>>b=[32;23;33;31];
>>[rank(B),rank([B,b])]
ans=
44
>>x=B\b
x=-5.4761
11.4934
-1.4292
2.4838
求
假设X=
>>x=B\b;x1=[1;1;1;1];X=x-x1
X=
-6.4761
10.4934
-2.4292
1.4838
求
>>norm(X)
ans=
12.6552
12.655就是
。
%求
>>norm(X)/norm(x1)
ans=
6.3276
6.3276即为
>>norm(a)
ans=
0.2244
>>norm(A)
ans=
30.2887
>>norm(a)/norm(A)%求
ans=
0.0074
所以
=0.0074
inv(A)%求A的逆矩阵,下求
>>d=inv(A);
>>norm(d)*norm(a)*norm(x)
ans=
288.4858
>>1-norm(d*(a))
ans=
-12.3693
>>288.4858/-12.3693
ans=
-23.3227
所以
=-23.322
>>norm(X)
ans=
0.0074
所以
>
(3)改变
,取a2=[00.0020.0010.003;0.0010.00400.001;0.003-0.004-0.0010;-0.001-0.0020-0.003]
B2=a2+A;C2=[Bb]
b=[32;23;33;31]
>>rank(B2)
ans=
4
>>rank(C2)
ans=
4
rank(B2)=rank(C2)
所以扰动矩阵有唯一解
>>x2=B2\b
x2=
1.0649
0.8893
1.0272
0.9859
>>x=B\b;x1=[1;1;1;1];X=x-x1%求
(设X=
)
X2=
-6.4761
10.4934
-2.4292
1.4838
norm(X2)%求
>>norm(X2)
ans=
12.6552
12.6552就是
>>norm(X)/norm(x1)%求
>>norm(X2)/norm(x1)
ans=
6.3276
所以
=6.3276
>>norm(a2)
ans=
0.0071
>>norm(a)/norm(A)%求
>>norm(a2)/norm(A)
ans=
2.3336e-004
所以
=2.3336e-004
norm(d)*norm(a2)*norm(x2)
ans=
9.0875
>1-norm(d*(a2))
ans=
0.6943
norm(X2)
ans=
12.6552
9.0875/0.6943
ans=
13.0887
所以
=13.0887
<
三、解三对角线性方程组的追赶法及其应用
⑴编写解三对角线性方程组的追赶法的通用程序,并应用于方程组
,检验程序的正确性;(解为
)
⑵求微分方程边值问题
的数值解(取步长
),并与精确解比较(精确解为
)。
说明:
离散化微分方程时,
解:
clearall;
a=[2,2,2,2,2];
b=[-1,-1,-1,-1];
c=[-1,-1,-1,-1];
r=[1,0,0,0,0];
n=length(a);
b=[0,b];
u
(1)=r
(1)/a
(1);
v
(1)=c
(1)/a
(1);
fork=2:
n-1
u(k)=(r(k)-u(k-1)*b(k))/(a(k)-v(k-1)*b(k));
v(k)=c(k)/(a(k)-v(k-1)*b(k));
end
u(n)=(r(n)-u(n-1)*b(n))/(a(n)-v(n-1)*b(n));
x(n)=u(n);
fork=n-1:
-1:
1
x(k)=u(k)-v(k)*x(k+1);
end
fprintf('Èý¶Ô½Ç·½³Ì×éµÄ½âΪ\n')
fork=1:
nfprintf('x(%1d)=%10.8f\n',k,x(k))
end
>>zhuiganfa%调用追赶法
三对角方程组的解为
x
(1)=0.83333333
x
(2)=0.66666667
x(3)=0.50000000
x(4)=0.33333333
x(5)=0.16666667
和结果
很好地吻合。
4、公元1225年,比萨的数学家Leonardo研究了方程
,得到一个根
,没有人知道他用什么方法得到这个值。
对于这个方程,分别用下列方法:
⑴迭代法
;⑵迭代法
;⑶对⑴的Steffensen加速方法;⑷对⑵的Steffensen加速方法;⑸Newton法。
求方程的根(可取
),计算到Leonardo所得到的准确度。
解:
由题意编写m文件如下
function[x0,k,er,x]=diedai(g,x0,wucha,max)
%g是给定的迭代函数
%x0是给定的初始值
%wucha是规定的误差范围
%max是所应许的最大迭代次数
%k是迭代次数加1
%x是不动点近似值
%x(x1,x2...,xn)
X
(1)=x0;
fork=2:
max
X(k)=feval('g',X(k-1));
k,er=abs(X(k)-X(k-1))
x=X(k);
if(er break;end ifk==max disp('超出迭代次数'); end end 其中定义的g函数是 functiony=g(x);y=20/(x^2+2*x+10); 在命令窗口中输入>>diedai('g',1,10^(-9),15) k= 15 er= 1.410125245193683e-005 超出迭代次数 X= Columns1through3 1.000000000000001.538461538461541.29501915708812 Columns4through6 1.401825309448601.354209390404291.37529809248738 Columns7through9 1.365929788170651.370086003401821.36824102361284 Columns10through12 1.369059812007481.368696397555521.36885768862873 Columns13through15 1.368786102577991.368817874396091.36880377314363 ans=1 取解为1.36880377314363 对于 ,只需修改diedai.m文件中的g,把其改为g1,编写m文件functiony=g1(x); y=(20-2*x^2-x^3)/10; 在命令窗口中输入diedai('g1',1,10^(-9),30) ….. k= 24 er= 1.36601568861328 k= 25 er= 1.36860974051282 k= 26 er= 1.36942327571766 取解为1.36860974051282 牛顿法: 编写m文件 function[p1,er,k,y]=ndf(f,df,p0,tol,max) %f是非线性函数 %df是f的微商 %p0是初始值 %tol是给定的允许误差 %max是迭代的最大次数 %p1是牛顿法求得的近似解 %er是p0的误差估计 %k是迭代次数 %y=f(p1) p0,feval('f',p0) fork=1: max p1=p0-feval('f',p0)/feval('df',p0); er=abs(p1-p0); p0=p1; p1,er,k,y=feval('f',p1) if((er break,end end 定义函数m文件 functiony=f(x) y=x^3+2*x^2+10*x-20; 定义微商函数m文件 functiony=df(x) y=3*x^2+4*x+10; 在命令窗口输入 >>ndf('f','df',1,10^(-9),10) p0=k= 11 ans=y=0.91756564217382 -7 p1=p1=1.411764705882351.36933647058824 er=er= 0.411764705882350.04242823529412 k=k= 24 y=y= 0.011148119412453.907985046680551e-014 p1=p1= 1.368808188617531.36880810782137 K=er= 31.776356839400251e-015 y=k= 1.704487072373695e-0065 p1=y= 1.368808107821370 er=ans= .0796******** 由结果知道牛顿法迭代到第三次已经达到所要求的精度 故方程的根为1.36880810782137 3. 编写Steffensen.m文件 i=2;x0=1;%设初始值 f=inline('20/(x^2+2*x+10)');%迭代格式 y=f(x0); z=f(y); x1=x0-(y-x0).^2/(z-2*y+x0); S.result=[x0;x1]; whileabs(x1-x0)>=1e-9%迭代精度 x0=x1; y=f(x0); z=f(y); x1=x0-(y-x0).^2/(z-2*y+x0); i=i+1; S.result(i)=x1; end S.step=[(0: i-1)]'; fprintf('迭代步数为: \t%d\n',i-1); forj=1: i fprintf('%10d',S.step(j));fprintf('\t'); f 在命令窗口输入Steffensen 迭代步数为: 4 01.0000000 11.3708139 21.3688082 31.3688081 41.3688081 分析结果知,Steffensen迭代加速步数减少了,到第四步已经达到了精度要求。 4.把迭代格式改为(20-2*x^2-x^3)/10,保存,在命令窗口输入 Steffensen 迭代步数为: 5 01.0000000 11.3334921 21.3684154 31.3688081 41.3688081 51.3688081 分析结果知,Steffensen迭代加速步数减少了,到第五步已经达到了精度要求。 五、用不同的数值方法计算积分 的近似值,其中 ⑴取不同的步长 ,分别用复合梯形公式和复合辛普森公式计算积分,比较两个公式的计算效果,是否存在一个最小的 ,使得精度不能再被改善? ⑵用龙贝格求积公式,取 ,并打印出T-表。 解: (1)用复合梯形公式,编写fhtx.m文件 functions=fhtx(f,a,b,n) %f是被积函数 %ab分别为积分的上下限 %n是子区间的个数 %s是梯形总面积 h=(b-a)/n; s=0; fork=1: (n-1) x=a+k*h; s=s+feval('f',x); end s=h*(feval('f',a)+feval('f',b))/2+h*s; 编写被积函数文件tf.m functiony=f(x); y=(10/x)^2*sin(10/x); 在命令窗口输入 >>fhtx('tf',1,3,10) ans= -4.7789 >>fhtx('tf',1,3,15) ans= -2.8604 >>fhtx('tf',1,3,20) ans= -2.2208 用复化辛普森公式,编写fhxps.m文件 functions=fhxps(f,a,b,n) %f是被积函数 %ab分别为积分的上下限 %n是子区间的个数 %s是梯形总面积 h=(b-a)/(2*n); s1=0;s2=0; fork=1: n x=a+(2*k-1)*h; s1=s1+feval('f',x); end fork=1: n-1 x=a+2*k*h; s2=s2+feval('f',x); end s=h*(feval('f',a)+feval('f',b)+4*s1+2*s2)/3; 在命令窗口输入 >>fhxps('f',1,3,15) ans= -1.4136 >>fhxps('f',1,3,20) ans= -1.4220 取n较大时 >>fhxps('f',1,3,1000) ans= -1.42602475569236 >>fhtx('tf',1,3,1000) ans= -1.42633620719670 >>fhtx('tf',1,3,2000) ans= -1.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 方法 计算 实习 题信计 101 和数