信计08级数值方法计算实习题2讲解Word下载.docx
《信计08级数值方法计算实习题2讲解Word下载.docx》由会员分享，可在线阅读，更多相关《信计08级数值方法计算实习题2讲解Word下载.docx（16页珍藏版）》请在冰豆网上搜索。
2.15
2.05
2.25
7.0
8.0
9.2
10.5
11.3
11.6
12
12.6
13.0
13.3
2.3
1.95
1.4
0.7
0.6
0.5
0.4
0.25
编写三次样条函数m文件：
function[f,f0]=scyt（x,y,y2_1,y2_N,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=followup（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）;
其中的followup函数为追赶法，其程序为
functionx=followup（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命令窗口输入：
>
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];
[f,f0]=scyt（x,y,0,0,x0,3.5）%其中给出了x=3.5处y的值
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
三次样条插值图像为>
plot（x,y,'
r*'
grid
拉格朗日插值函数程序，定义m文件
functionf=lglr（x,y,x0）
l=y（i）;
for（j=1:
i-1）;
l=l*（t-x（j））/（x（i）-x（j））;
end;
f=f+l;
simplify（f）;
if（i==n）
if（nargin==3）
f=subs（f,'
x0）
else
f=collect（f）;
f=vpa（f,6）;
%插值多项式的系数化为6位小数
end
在命令窗口输入>
f=lglr（x,y）
-337.863*t^2+219.386*t-.171512e-3*t^13+.706734*t^8+.152386e-4*t^14-.105677e-5*t^15-.271334*t^9+304.616*t^3-.108815e-1*t^11+.622166e-1*t^10+65.2964*t^5+.152863e-2*t^12-13.7172*t^6+.477508*t^7+.559234e-7*t^16+.769797e-14*t^20-.985682e-12*t^19-60.7814-176.003*t^4+.588899e-10*t^18-.217905e-8*t^17%拉格朗日20次公式
f=lglr（x,y,3.5）
194.2495
二、已知Wilson矩阵，且向量，则方程组有准确解。
⑴用Matlab内部函数求，的所有特征值和；
:
A=[10787;
7565;
86109;
75910];
B=det（A）
B=
1
A的所以特征值：
D=eig（A）
D=
0.0102
0.8431
3.8581
30.2887
条件数：
30.2887/0.0102
ans=
2.9695e+003
Cond（A）=Y1/Yn=30.2887/0.0102=2.9695e+003
⑵令，解方程组，并求出向量和，从理论结果和实际计算结果两方面分析方程组解的相对误差与的相对误差的关系；
⑶再改变扰动矩阵（其元素的绝对值不超过0.005），重复第2问。
三、解三对角线性方程组的追赶法及其应用
⑴编写解三对角线性方程组的追赶法的通用程序，并应用于方程组
，检验程序的正确性；
（解为）
解：
由最赶法的递推公式（课本p.160）
编写m文件
a=[2,2,2,2,2];
b=[-1,-1,-1,-1];
c=[-1,-1,-1,-1];
f=[1,0,0,0,0];
n=length（a）;
b=[0,b];
u
（1）=f
（1）/a
（1）;
v
（1）=c
（1）/a
（1）;
fork=2:
u（k）=（f（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:
formatrat
fprintf（'
x（%1d）=%10.8f\n'
k,x（k））
end
并且以zgf保存m文件
在命令窗口中输入>
zgf
三对角方程组的解为
x
（1）=0.83333333
x
（2）=0.66666667
x（3）=0.50000000
x（4）=0.33333333
x（5）=0.16666667
和结果很好地吻合。
⑵求微分方程边值问题的数值解（取步长），并与精确解比较（精确解为）。
说明：
离散化微分方程时，
四、公元1225年，比萨的数学家Leonardo研究了方程，得到一个根，没有人知道他用什么方法得到这个值。
对于这个方程，分别用下列方法：
⑴迭代法；
⑵迭代法；
⑶对⑴的Steffensen加速方法；
⑷对⑵的Steffensen加速方法；
⑸Newton法。
求方程的根（可取），计算到Leonardo所得到的准确度。
1解：
由题意编写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
超出迭代次数'
X
其中定义的g函数是
functiony=g（x）;
y=20/（x^2+2*x+10）;
diedai（'
1,10^（-9）,15）
...
k=
15
er=
1.410125245193683e-005
超出迭代次数
X=
Columns1through3
1.00000000000000
升级会员 