计算方法.docx
《计算方法.docx》由会员分享，可在线阅读，更多相关《计算方法.docx（29页珍藏版）》请在冰豆网上搜索。
计算方法
《计算方法》实验指导书
二级学院：
专业：
指导教师：
班级学号：
姓名：
实验一插值与拟合
1、实验目的：
①通过用Matlab编程和插值与拟合中的某种具体算法解决具体问题，更深一步的体会数值分析这门课的重要性，同时加深对插值与拟合公式某种具体算法的理解。
②熟悉Matlab编程环境。
2、实验要求:
用Matlab和插值与拟合中的某种具体算法编写代码并执行，完成解决具体问题。
3、实验内容：
4、题目：
5、原理：
6、设计思想：
7、对应程序：
8、实验结果：
9、图形（如果可视化）
10、实验体会：
11、样例
实验一牛顿插值
实验目的：
①通过用Matlab编程解决数值分析问题，更深一步的体会数值分析这门课的重要性，同时加深对牛顿插值公式的理解。
②熟悉Matlab编程环境。
实验要求:
用Matlab编写代码并执行，完成第一章34题（修改）。
实验内容：
题目：
构造适合下列数据表的牛顿插值公式：
x
-1
0
1
3
y
-1
1
3
5
6
1
原理：
牛顿三次插值公式P3=f（x0）+f[x0,x1]（x-x0）+f[x0,x1,x2]（x-x0）（x-x1）+f[x0,x1,x2,x3]（x-x0）（x-x1）（x-x2）
设计思想：
由于牛顿三次插值公式是由多项式组成的，且存在极大的规律性：
首先，从第二项看起，都是差商与x的表达式的乘积，所以有：
P3=Y
（1）;
a=0;
b=（x-X
（1））;
然后：
每循环一次，更新三者的值：
a=CS（i,1）;//CS为差商矩阵。
第一行为一阶差商，第二行为二阶差商，第三行为三阶差商。
（左三角阵）.
P3=P3+a*（b）;
b=（b）*（x-X（i+1））;
源代码：
（图）
对应程序：
（建议用源文件来运行，因为直接粘贴，存在多文件和文件名的问题）
disp（'第一章34题：
构造适合下列数据表的牛顿插值公式：
'）
disp（'x-1013'）
disp（'y-1135'）
disp（'diff（y）61'）
%三次牛顿插值公式为：
%P3=f（x0）+f[x0,x1]（x-x0）+f[x0,x1,x2]（x-x0）（x-x1）+f[x0,x1,x2,x3]（x-x0）（x-x1）（x-x2）
%=1+2*x-1/12*（x+1）*x*（x-1）
X=[-1013];
Y=[-1135];
CS=[100,100,100;100,100,100;100,100,100];
fori=1:
1:
3
CS（1,i）=（Y（i+1）-Y（i））/（X（i+1）-X（i））;
end
fori=1:
1:
2
CS（2,i）=（CS（1,i+1）-CS（1,i））/（X（i+2）-X（i））;
end
CS（3,1）=（CS（2,2）-CS（2,1））/（X（4）-X
（1））;
CS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
symsP3
symsx
symsb
P3=Y
（1）;
a=0;
b=（x-X
（1））;
fori=1:
1:
3
a=CS（i,1）;
P3=P3+a*（b）;
b=（b）*（x-X（i+1））;
end
fprintf（'\n求得牛顿插值公式为：
'）
P3
fprintf（'/****************************************/'）
实验结果：
图形（牛顿插值公式）
实验体会：
本实验是第一次接触到这种新的软件。
在编程过程中，考虑了以前没有考虑过的问题，比如：
误差的传递。
感受到了选择方法的重要性。
实验一插值与拟合
1、实验目的：
①通过用Matlab编程和插值与拟合中的某种具体算法解决具体问题，更深一步的体会数值分析这门课的重要性，同时加深对插值与拟合公式某种具体算法的理解。
②熟悉Matlab编程环境。
2、实验要求:
用Matlab和插值与拟合中的某种具体算法编写代码并执行，完成解决具体问题。
3、实验内容：
4、题目：
5、原理：
6、设计思想：
7、对应程序：
8、实验结果：
9、图形（如果可视化）
10、实验体会：
11、样例
实验二辛甫生公式和梯形公式
实验目的：
①熟练编程技巧，体会辛甫生公式和梯形公式的用法。
②理解Matlab下一种处理积分的方法。
实验要求:
编写程序，求出第二章17题（修改）。
要求给出实验结果。
实验内容：
题目：
已知y=
，用辛甫生公式和梯形公式分别求值PI=
.
原理：
1.辛甫生公式：
=
。
2．梯形公式：
=
。
设计思想：
首先，找到表达式中的共同点----一个函数。
编写该函数。
然后，分别编写两公式的表达式。
源代码：
（下图为程序的源代码）
对应程序：
（建议用源文件来运行，因为直接粘贴，存在多文件和文件名的问题）
disp（'第二章17题：
已知y=4/（1+x.^2）,'）
disp（'用辛甫生公式和梯形公式分别求积分值PI=int（y,x,0,1）'）
%辛甫生公式f（x）在[a,b]上的积分为：
（b-a）/6*（f（a）+f（b）+4*f（（a+b）/2））;
%梯形公式f（x）在[a,b]上的积分为：
（b-a）/2*（f（a）+f（b））;
fprintf（'\n方法一：
辛甫生公式计算\n'）
b=1;
a=0;
PI=（b-a）/6*（f（a）+f（b）+4*f（（a+b）/2））
fprintf（'\n方法二：
梯形公式计算\n'）
PI=（b-a）/2*（f（a）+f（b））
（下图为函数的源代码：
）
对应程序：
（保存为文件名f.m文件）
functiony=f（x）
y=4/（1+x.^2）;
实验结果：
实验体会：
单从计算结果就可以得出这样的结论：
梯形公式没有辛甫生公式的精度高。
它的相对误差大一些。
还一点就是，除了用差商的方法外，又学会了一种求积分的方法。
实验三
实验目的：
深入理解龙格—库塔方法的原理，同时学会用迭代方法看待某些特定问题。
比较本方法与其它方法（尤其是欧拉方法）解题的不同处，分析两者的精度。
实验要求:
编写程序，用龙格—库塔方法求解第三章3题（修改）
实验内容：
题目：
取h=0.1，用四阶龙格-库塔方法求解初值问题：
=
y（0）=0
并与精确解y=
比较结果.
原理：
本题很明显为迭代方法的一种。
已知初值和函数的导数。
再根据四阶龙格—库塔公式：
k1=dy（x0,y0）;
k2=dy（（2*x0+h）/2,y0+h/2*k1）;
k3=dy（（2*x0+h）/2,y0+h/2*k2）;
k4=dy（x0+h,y0+h*k3）;
y0=y0+h/6*（k1+2*k2+2*k3+k4）;
x0=x0+h;
设计思想：
由四阶龙格—库塔公式很容易看出，只要作一个循环即可完与此项操作。
每次更新k1,k2,k3,k4,y0以及x0的值。
这样即可得出y。
源代码：
（下图为程序的源代码）
对应程序：
（建议用源文件来运行，因为直接粘贴，存在多文件和文件名的问题）
disp（'第三章3题：
取h=0.1，用四阶龙格-库塔方法求解初值问题：
'）
disp（'diff（y,x,1）=1/（1+x.^2）-2*y.^21<=x<=4'）
disp（'y（0）=0'）
disp（'并与精确解y=x/（1+x.^2）比较结果.'）
x0=1;%初值为1
y0=0.5;%初值为0.5
h=input（'请输入步长>>'）
fori=1:
1:
5
k1=dy（x0,y0）;
k2=dy（（2*x0+h）/2,y0+h/2*k1）;
k3=dy（（2*x0+h）/2,y0+h/2*k2）;
k4=dy（x0+h,y0+h*k3）;
y0=y0+h/6*（k1+2*k2+2*k3+k4）;
x0=x0+h;
fprintf（'\nk%d=%f,k%d=%f,k%d=%f,k%d=%f',1,k1,2,k2,3,k3,4,k4）
fprintf（'\n=>>y[%d]=%f',i,y0）;
end
fprintf（'\n以下为用公式：
y=x/（1+x.^2）求得的精确值\n'）
x0=1;
fori=0:
1:
5
y0=f（x0）;
x0=x0+h;
fprintf（'y[%d]=%f',i,y0）;
end
fprintf（'\n/****************************************************/'）
（下图为yy=dy（x）以及y=f（x）两函数的源代码：
）
文件名：
f.m
functiony=f（x）
y=x./（1+x.^2）;
文件名：
dy.m
functionyy=dy（x,y）
yy=1/（1+x.^2）-2*y.^2;
实验结果：
图形
实验体会：
本实验中，关于龙格—库塔方法的使用，使我个人又对迭代方法有了新的一种理念。
相对于其它迭代方法而言，龙格库塔方法具有更高的精度。
但其缺点则是给人带来更大的计算量。
实验四
实验目的：
练习动手推算迭代分式的能力，加深理解迭代的真谛。
体会迭代与极限间的微妙联系。
实验要求:
编写一段代码，实验第四章6题。
实验内容：
题目：
给出计算x=
的迭代公式，讨论迭代过程的收敛性，并证明x=2
原理：
根据上式的规律可以看出，若将每对小括号内的表达式看作一个整体，则这个整体总会呈现出迭代的规律。
设计思想：
经以上分析，令x=（2+x.^0.5）作为迭代公式。
源代码：
（下图为程序的源代码）
对应程序：
（建议用源文件来运行，因为直接粘贴，存在多文件和文件名的问题）
disp（'第四章6题：
给出计算x=（2+（2+（2+（（2+...）^0.5）^0.5）^0.5）^0.5）^0.5'）
disp（'的迭代公式，讨论迭代过程的收敛性，并证明x=2'）
fprintf（'\n迭代公式为：
x=（2+x）^0.5\n'）
x1=2^0.5;
x2=2;
fazi=input（'请输入精度（当精度为0时，说明为精确值！
）>>'）
whileabs（x2-x1）>fazi
x1=（2+x1）^0.5;
x2=（2+x1）^0.5;
end
ifx1==2
fprintf（'\n此时为准确值：
x=%f\n',x1）
else
fprintf（'\n近似值：
x=%f',x1）
end
fprintf（'\n/*******************************************/'）
实验结果：
实验体会：
迭代的过程即刷新变量的过程。
在每次更新值动作中，变量会呈现出相当的规律性。
对于本题，还有一个很奇妙的地方，就是当选定精度为0，即没有误差时，恰好为该式的极限。
也是该式的精确解。
实验五
实验目的：
练习引入迭代和矩阵的形式来解决方程组问题。
巩固方程与矩阵相互关系的概念。
实验要求:
学会用矩阵的形式来解决方程组问题在计算机上的实现。
分析雅可比迭代与其它迭代的异同。
实验内容：
题目：
用雅可比迭代求方程组：
AX=B
已知：
A、B如下：
A=B=
-4-10-1000
-14-10-105
0-1400-10
-1004-106
0-10-14-12
00-10-146
原理：
先找出下三角阵-L，再找出上对角阵-U，还有主对角阵D。
迭代公式x=inv（D）*（L+U）*x+inv（D）*B
设计思想：
利用迭代即可在循环中实现的原则来完成此方程组的求解，所用的迭代公式为：
x=inv（D）*（L+U）*x+inv（D）*B
源代码：
（下图为程序的源代码）
对应程序：
（建议用源文件来运行，因为直接粘贴，存在多文件和文件名的问题）
disp（'第五章4题：
用雅可比迭代求方程组：
AX=B'）
fprintf（'已知：
A、B如下：
'）
A=[-4-10-100;-14-10-10;0-1400-1;-1004-10;0-10-14-1;00-10-14]
X=[0;0;0;0;0;0];%X初值为[0;0;0;0;0;0]
B=[0;5;0;6;-2;6]
D=[-400000;040000;004000;000400;000040;000004];
L=-[000000;-100000;0-10000;-100000;0-10-100;00-10-10];
U=-[0-10-100;00-10-10;00000-1;0000-10;00000-1;000000];
fori=1:
1:
5
X=inv（D）*（L+U）*X+inv（D）*B;
end
fprintf（'\n方程组的解为：
'）
X
fprintf（'\n/**************************************************/'）
实验结果：
实验体会：
利用Matlab编写程序计算方程组的问题是很容易的，但也存在着一定的问题，那就是精度问题，在系数相差很大时，方程组解往往出现很大误差。
实验六
实验目的：
加强编程能力和编程技巧，练习从数值分析的角度看问题。
同时用Matlab编写代码。
实验要求:
用选主元素法和高斯消去法两种方法解方程组。
实验内容：
题目：
第六章3题（修改）
>>第六章3题：
用选主元素法和高斯消去法求解下列方程组：
已知方程增广矩阵：
a=
2-1006
-1-3-201
-13-200
00-351
原理：
（1）高斯消去法：
相对于约当消去法而言，总的来说就是将增广矩阵化为下三角阵。
（2）选主元素法:
相对于高斯消去法的唯一不同是每次将对角线元素化为1前，要先按当前要排列元素所在列大小进行排列。
设计思想：
（1）高斯消去法：
首先让每一行的元素除以该行的主对角线元素。
然后利用此行使位于下一行主对角线以前的元素变为0，依次类推。
（2）选主元素法：
在高斯消去法的基础上，每次进行化上三角阵之前，重新排列各方程的位置。
源代码：
（下图为高斯消去法程序的源代码）
对应程序：
（建议用源文件来运行，因为直接粘贴，存在多文件和文件名的问题）
（1）高斯消去法
disp（'第六章3题：
用高斯消去法求解下列方程组：
'）
a=[2-1006;-1-3-201;-13-200;00-351]
i=1;
z=2;
forn=[321]
ifn==3%y用来控制下面循环的次数
y=3;
elseifn==2
y=2;
elseifn==1
y=1;
end%y用来控制下面循环的次数
h=z;
forx=（1:
y）%%此段用来消元
d=diag（a）;
c=a（:
i）;
ifc（h）~=0%用来判断是否为0
a（h,:
）=（-d（i）/c（h）*a（h,:
））+a（i,:
）;
end%用来判断是否为0
h=h+1;
end%%此段用来消元
z=z+1;
i=i+1;
d=diag（a）;
end
forn=（1:
4）%把主对角的元素全部化为1
a（n,:
）=a（n,:
）/d（n）;
end%把主对角的元素全部化为1
x=3;%此段把得到的矩阵化为行最简形矩阵
forn=（1:
3）
y=4;
fork=（1:
n）
u=a（x,:
）;
p=u（y）;
a（x,:
）=a（x,:
）-p*a（y,:
）;
y=y-1;
end%此段把得到的矩阵化为行最简形矩阵
x=x-1;
end
fprintf（'最后得到的矩阵是：
'）
a
fprintf（'用高斯法解方程得到的解是\n'）
b=a（:
5）
fprintf（'/*****************************************************/'）
（下图为选主元素法程序的源代码）
对应程序：
（建议用源文件来运行，因为直接粘贴，存在多文件和文件名的问题）
（选主元素法）
disp（'第六章3题：
用主元消去法求解下列方程组：
'）
a=[2-1006;-1-3-201;-13-200;00-351]
i=1;
z=2;
forn=[321]
c=a（:
i）;%选主元素
d=diag（a）;
j=1;
forww=（1:
n）
ifabs（d（i））yy=a（i,:
）;
a（i,:
）=a（n+1-j,:
）;
a（n+1-j,:
）=yy;
end
j=j+1;
end%选主元素
ifn==3%y用来控制下面循环的次数
y=3;
elseifn==2
y=2;
elseifn==1
y=1;
end%y用来控制下面循环的次数
h=z;%%此段用来消元
forx=（1:
y）
d=diag（a）;
c=a（:
i）;
ifc（h）~=0%用来判断是否为0
a（h,:
）=（-d（i）/c（h）*a（h,:
））+a（i,:
）;
end%用来判断是否为0
h=h+1;
end%%此段用来消元
z=z+1;
i=i+1;
d=diag（a）;
end
forn=（1:
4）%把主对角的元素全部化为1
a（n,:
）=a（n,:
）/d（n）;
end%把主对角的元素全部化为1
x=3;%此段把得到的矩阵化为行最简形矩阵
forn=（1:
3）
y=4;
fork=（1:
n）
u=a（x,:
）;
p=u（y）;
a（x,:
）=a（x,:
）-p*a（y,:
）;
y=y-1;
end
a;
x=x-1;%此段把得到的矩阵化为行最简形矩阵
end
fprintf（'最后得到的矩阵是：
'）
a
fprintf（'用选主元法解方程得到的解是:
'）
b=a（:
5）
fprintf（'/********************************************************/'）
实验结果：
实验体会：
主元消去法和高斯消去法的确是两个非常锻炼人编程的方法。
在这两个程序中，循环是必不可少的，也就是说要更新大量的变量。
不过，主元消去法也好，高斯消去法也好，它们的原理还算不难理解。
执笔人：
任玉杰审核人：
于华
升级会员 