数值方法计算实习题1文档格式.docx
《数值方法计算实习题1文档格式.docx》由会员分享,可在线阅读,更多相关《数值方法计算实习题1文档格式.docx(23页珍藏版)》请在冰豆网上搜索。
%插值点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
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);
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)
Error:
对角有元素为0!
'
end;
d=ones(n,1);
a=ones(n-1,1);
c=ones(n-1);
n-1)
a(i,1)=A(i+1,i);
c(i,1)=A(i,i+1);
d(i,1)=A(i,i);
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);
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);
在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)
else
l=y(i);
for(j=1:
i-1)
l=l*(t-x(j))/(x(i)-x(j));
end;
for(j=i+1:
f=f+l;
simplify(f);
if(i==n)
if(nargin==3)
f=subs(f,'
f=collect(f);
f=vpa(f,6);
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];
>
[rank(B),rank([B,b])]
44
x=B\b
x=-5.4761
11.4934
-1.4292
2.4838
求
假设X=
x=B\b;
x1=[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)
6.3276
6.3276即为
norm(a)
0.2244
norm(A)
norm(a)/norm(A)%求
0.0074
所以
=0.0074
inv(A)%求A的逆矩阵,下求
d=inv(A);
norm(d)*norm(a)*norm(x)
288.4858
1-norm(d*(a))
-12.3693
288.4858/-12.3693
-23.3227
=-23.322
(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]
31]
rank(B2)
4
rank(C2)
rank(B2)=rank(C2)
所以扰动矩阵有唯一解
x2=B2\b
x2=
1.0649
0.8893
1.0272
0.9859
X=x-x1%求
(设X=
X2=
1.4838
norm(X2)%求
norm(X2)
12.6552就是
norm(X)/norm(x1)%求
norm(X2)/norm(x1)
=6.3276
norm(a2)
0.0071
norm(a2)/norm(A)
2.3336e-004
=2.3336e-004
norm(d)*norm(a2)*norm(x2)
9.0875
1-norm(d*(a2))
0.6943
norm(X2)
9.0875/0.6943
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:
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));
u(n)=(r(n)-u(n-1)*b(n))/(a(n)-v(n-1)*b(n));
x(n)=u(n);
fork=n-1:
1
x(k)=u(k)-v(k)*x(k+1);
fprintf('
È
ý
¶
Ô
½
Ç
·
³
Ì
×
é
µ
Ä
â
Î
ª
\n'
fork=1:
nfprintf('
x(%1d)=%10.8f\n'
k,x(k))
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;
max
X(k)=feval('
g'
X(k-1));
k,er=abs(X(k)-X(k-1))
x=X(k);
if(er<
wucha),
ifk==max
超出迭代次数'
end
end
其中定义的g函数是
functiony=g(x);
y=20/(x^2+2*x+10);
在命令窗口中输入>
diedai('
1,10^(-9),15)
k=
15
er=
1.410125245193683e-005
超出迭代次数
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)
…..
24
1.36601568861328
25
1.36860974051282
26
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)
p1=p0-feval('
p0)/feval('
df'
p0);
er=abs(p1-p0);
p0=p1;
p1,er,k,y=feval('
p1)
if((er<
tol)||(y==0)),
break,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('
'
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;
S.step=[(0:
i-1)]'
;
迭代步数为:
\t%d\n'
i-1);
forj=1:
i
fprintf('
%10d'
S.step(j));
\t'
f
在命令窗口输入Steffensen
01.0000000
11.3708139
21.3688082
31.3688081
41.3688081
分析结果知,Steffensen迭代加速步数减少了,到第四步已经达到了精度要求。
4.把迭代格式改为(20-2*x^2-x^3)/10,保存,在命令窗口输入
Steffensen
5
11.3334921
21.3684154
51.3688081
分析结果知,Steffensen迭代加速步数减少了,到第五步已经达到了精度要求。
五、用不同的数值方法计算积分
的近似值,其中
⑴取不同的步长
,分别用复合梯形公式和复合辛普森公式计算积分,比较两个公式的计算效果,是否存在一个最小的
,使得精度不能再被改善?
⑵用龙贝格求积公式,取
,并打印出T-表。
(1)用复合梯形公式,编写fhtx.m文件
functions=fhtx(f,a,b,n)
%f是被积函数
%ab分别为积分的上下限
%n是子区间的个数
%s是梯形总面积
h=(b-a)/n;
s=0;
(n-1)
x=a+k*h;
s=s+feval('
x);
s=h*(feval('
a)+feval('
b))/2+h*s;
编写被积函数文件tf.m
functiony=f(x);
y=(10/x)^2*sin(10/x);
fhtx('
tf'
1,3,10)
-4.7789
1,3,15)
-2.8604
1,3,20)
-2.2208
用复化辛普森公式,编写fhxps.m文件
functions=fhxps(f,a,b,n)
h=(b-a)/(2*n);
s1=0;
s2=0;
x=a+(2*k-1)*h;
s1=s1+feval('
x=a+2*k*h;
s2=s2+feval('
b)+4*s1+2*s2)/3;
fhxps('
-1.4136
-1.4220
取n较大时
1,3,1000)
-1.42602475569236
-1.42633620719670
1,3,