计算实习MATLAB编程1.docx
《计算实习MATLAB编程1.docx》由会员分享，可在线阅读，更多相关《计算实习MATLAB编程1.docx（48页珍藏版）》请在冰豆网上搜索。
计算实习MATLAB编程1
计算实习（MATLAB编程）
第一章
1.用3位数字计算出方程：
的解x，y，再用6位数字计算出x与y，已知正确解为练习练习x=1,y=-1,计算结果说明什么？
[x,y]=solve（'0.780*x+0.563*y=0.217','0.457*x+0.330*y=0.127'）
x=
1.
y=
-1.
a=[0.7800.563;0.4570.330];b=[0.217;0.127];
l=vpa（a（2,1）/a（1,1）,3）;
c=vpa（a（2,2）-a（1,2）*l,3）;d=vpa（b
（2）-b
（1）*l,3）;
y=vpa（d/c,3）,x=vpa（（b
（1）-a（1,2）*y）/a（1,1）,3）
y=
-1.98
x=
1.71
a=[0.7800.563;0.4570.330];b=[0.217;0.127];
l=vpa（a（2,1）/a（1,1）,6）;
c=vpa（a（2,2）-a（1,2）*l,6）;d=vpa（b
（2）-b
（1）*l,6）;
y=vpa（d/c,6）,x=vpa（（b
（1）-a（1,2）*y）/a（1,1）,6）
y=
-.997571
x=
.998247
第二章
用MATLAB作插值计算：
设n个节点以数组x0,y0输入，m个插值点以数组x输入，输出数组y为m个插值。
比如一个名为lagr1.m的M文件内容为：
functiony=lagr1（x0,y0,x）
n=length（x0）;m=length（x）;
fori=1:
m
z=x（i）;
s=0.0;
for k=1:
n
p=1.0;
forj=1:
n
if j~=k
p=p*（z-x0（j））/（x0（k）-x0（j））;
end
end
s=s+p*y0（k）;
end
y（i）=s;
end
作Lagrange插值时只需在输入x0,y0,x后运行：
y=lagr1（x0,y0,x）即可；
分段线性插值有现成的程序：
y=interp1（x0,y0,x）
三次样条插值也有现成的程序：
y=interp1（x0,y0,x,'spline'）;此种方法属于自然边界条件。
其他边界条件下的三次样条插值程序见MATLAB的样条工具包（SplineToolbook）。
例：
对，，用个节点（等分）作上述三种插值，用m（=21）个插值点（等分）作图，比较结果。
其语句如下：
n=11;m=21;
x=-5:
10/（m-1）:
5;
y=1./（1+x.^2）;
z=0*x;
x0=-5:
10/（n-1）:
5;
y0=1./（1+x0.^2）;
y1=lagr1（x0,y0,x）;
y2=interp1（x0,y0,x）;
y3=interp1（x0,y0,x,'spline'）;
[x'y'y1'y2'y3']
plot（x,z,'r',x,y,'k:
',x,y1,'y',x,y2,'b',x,y3,'o'）
ans=
-5.0000 0.0385 0.0385 0.0385 0.0385
-4.5000 0.0471 1.5787 0.0486 0.0484
-4.0000 0.0588 0.0588 0.0588 0.0588
-3.5000 0.0755 -0.2262 0.0794 0.0745
-3.0000 0.1000 0.1000 0.1000 0.1000
-2.5000 0.1379 0.2538 0.1500 0.1401
-2.0000 0.2000 0.2000 0.2000 0.2000
-1.5000 0.3077 0.2353 0.3500 0.2973
-1.0000 0.5000 0.5000 0.5000 0.5000
-0.5000 0.8000 0.8434 0.7500 0.8205
0 1.0000 1.0000 1.0000 1.0000
0.5000 0.8000 0.8434 0.7500 0.8205
1.0000 0.5000 0.5000 0.5000 0.5000
1.5000 0.3077 0.2353 0.3500 0.2973
2.0000 0.2000 0.2000 0.2000 0.2000
2.5000 0.1379 0.2538 0.1500 0.1401
3.0000 0.1000 0.1000 0.1000 0.1000
3.5000 0.0755 -0.2262 0.0794 0.0745
4.0000 0.0588 0.0588 0.0588 0.0588
4.5000 0.0471 1.5787 0.0486 0.0484
5.0000 0.0385 0.0385 0.0385 0.0385
插值与拟合实验内容：
1． 选择一些函数，在n个节点上，用拉格朗日、分段线性、三次样条三种插值方法，计算m个插值点的函数值，通过数值和图形输出。
a.y=sinx,
n=10;m=100;
x=linspace（0,2*pi,100）;
y=sin（x）;
z=x*0;
x0=linspace（0,2*pi,6）;
y0=sin（x0）;
y1=lagr1（x0,y0,x）;
y2=interp1（x0,y0,x）;
y3=interp1（x0,y0,x,'spline'）;
plot（x,z,x,y,'k',x,y1,x,y2,x,y3,'y'）
b.
n=10;m=100;b=0:
0.1:
1;a=b*0;
x=linspace（-1,1,100）;
y=sqrt（1-x.*x）;
z=x*0;
x0=linspace（-1,1,6）;
y0=sqrt（1-x0.*x0）;
y1=lagr1（x0,y0,x）;
y2=interp1（x0,y0,x）;
y3=interp1（x0,y0,x,'spline'）;
plot（x,z,x,y,'k',x,y1,x,y2,x,y3,'b',a,b,'k'）;
2． 用在x=0,1,4,9,16产生5个节点，用不同的节点构造插值公式来计算x=5处的插值，与精确值比较并进行分析。
x0=[014 916];y0=sqrt（x0）;x=5;y=sqrt（5）
y1=lagr1（x0,y0,x）
x0
（1）=[];y0=sqrt（x0）;y2=lagr1（x0,y0,x）
x0（4）=[];y0=sqrt（x0）;y3=lagr1（x0,y0,x）
x0
（1）=[];y0=sqrt（x0）;y4=lagr1（x0,y0,x）
y=
2.2361
y1=
2.0794
y2=
2.2540
y3=
2.2667
y4=
2.2000
第三章
1901年龙格（Runge）给出一个例子：
，定义在区间[-1,1]上，这是一个很光滑的函数，它的任意阶导数都存在，对它在[-1,1]上作等距节点插值时，插值多项式的情况见图1
从图1中看出，在靠近-1或1时余项会随着n的增大而很大，譬如，但，从图中还可看到在0附近插值效果是好的，即余项较小。
另一个现象是插值多项式随节点增多而振动更多。
这种插值多项式当节点增加时不能更好地接近被插值函数的现象称为龙格现象。
n=4;m=8;k=12;
x=-1:
0.01:
1;
y=1./（1+25*x.^2）;
z=0*x;
x0=-1:
2/n:
1;
y0=1./（1+25*x0.^2）;
z0=lagr1（x0,y0,x）;
x1=-1:
2/m:
1;
y1=1./（1+25*x1.^2）;
z1=lagr1（x1,y1,x）;
x2=-1:
2/k:
1;
y2=1./（1+25*x2.^2）;
z2=lagr1（x2,y2,x）;
plot（x,y,x,z0,x,z1,x,z2,x,z,'y'）
legend（'f（x）','n=4','n=8','n=12'）
第四章
例2 为了开拓市场，某公司对其新产品作了一系列调查，他们发现这一新产品的销量与下列事件关系密切：
其一是温度，其二是上证指数，其三是广告费……为此，他们记录了下面的数据表，假设这些事件与销售量近似成线性关系，试给出这种关系的数学表达式。
记录
温度
上证指数
广告费
推销员数
返修率
销售量
1
39
567
10000
2
0.20
75
2
37
679
0
3
0.15
68
3
30
346
5000
3
0.10
105
4
25
987
5000
3
0.08
136
5
25
1101
0
4
0.07
152
6
10
1004
5000
5
0.07
191
7
15
667
0
4
0.05
148
8
6
604
10000
6
0.04
234
解 对于这个例子，我们可以把销售量作为y，并用分别表示温度、上证指数、广告费、推销人员数目、返修率。
然后，利用最小二乘法就可以得到其关系表达式：
A=[39,567,10000,2,0.20;37,679,0,3,0.15;30,346,5000,3,0.10;25,987,5000,3,0.08;
25,1101,0,4,0.07;10,1004,5000,5,0.07;15,667,0,4,0.05;5,604,10000,6,0.04];
b=[75;68;105;136;152;191;148;234];
a=A\b
a=
0.6877
0.0437
0.0049
29.6072
-447.3589
a1=[0.147;0.051;4.685;28.28;-351.839];
norm（A*a-b）
norm（A*a1-b）
ans=
16.7645
ans=
7.7607e+004
第六章 快速傅立叶变换
1试用高级语言写出计算周期函数三角插值的子程序。
并分别对N=16，64，128计算周期三角函数的系数ak和bk，其中周期三角函数的波形图如图7所示。
x=inline（'pi-abs（t）'）;
t=-pi:
0.1:
pi;
T=x（t）;
y1=0;x1=0:
0.01:
pi;
plot（t,T,t,y1,0,x1）;
图7
其MATLAB程序为：
N=input（'N=?
'）；
x=inline（'
升级会员 