数值分析实验报告(Matlab实现).doc
《数值分析实验报告(Matlab实现).doc》由会员分享，可在线阅读，更多相关《数值分析实验报告(Matlab实现).doc（25页珍藏版）》请在冰点文库上搜索。
学生实验报告
实验课程名称数值分析
开课实验室数学与统计学院实验室
学院2010年级数学与应用数学专业班01班
学生姓名学号
开课时间2012至2013学年第一学期
总成绩
教师签名
课程
名称
数值分析
实验项目
名称
Gauss消元法
实验项目类型
验证
演示
综合
设计
其他
指导
教师
何光辉
成绩
√
一、实验目的：
（1）高斯列主元消去法求解线性方程组的过程
（2）熟悉用迭代法求解线性方程组的过程
（3）设计出相应的算法，编制相应的函数子程序
二、实验内容
分别用高斯列主元消元法和直接消元法求解线性方程组：
三、实验原理
对于线性方程组
（1）
常记为矩阵形式
（2）
根据高等代数的知识，若，上式的解存在且唯一。
（1）Gauss直接消元法
考虑上述线性方程组的增广矩阵，对增广矩阵进行行变换，将
（2）式化为等价的三角形方阵，然后回代解之，这就是Gauss消元法。
具体如下：
a）消元
①令；
②对k=1到n-1，若，进行
b）回代，若
（2）Gauss列主元消元法
设列主元消元法已完成的第k-1（）次消元，的到方程组
在进行第k次消元前，先进行2个步骤：
a）在至这一列内选出最大值，即，若，此时
方程组无确定解，应给出退出信息。
b）若，则交换第行和行，然后用Gauss消元法进行消元。
四、MATLAB软件实现
（1）写出Gauss消元法和列主元消元法实现的MATLAB函数
根据以上的算法，写出如下程序：
%%%%%%%Gauss消元法%%%%%%%%%%%%
functiony=Gauss1（A,b）
[m,n]=size（A）;
%检查系数正确性
ifm~=n
error（'矩阵A的行数和列数必须相同'）;
return;
end
ifm~=size（b）
error（'b的大小必须和A的行数或A的列数相同'）;
return;
end
%再检查方程是否存在唯一解
ifrank（A）~=rank（[A,b]）
error（'A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解'）;
return;
end
%这里采用增广矩阵行变换的方式求解
c=n+1;
A（:
c）=b;
%%消元过程
fork=1:
n-1
A（k+1:
n,k:
c）=A（k+1:
n,k:
c）-（A（k+1:
n,k）/A（k,k））*A（k,k:
c）;
end
%%回代结果
x=zeros（length（b）,1）;
x（n）=A（n,c）/A（n,n）;
fork=n-1:
-1:
1
x（k）=（A（k,c）-A（k,k+1:
n）*x（k+1:
n））/A（k,k）;
end
%显示计算结果
%disp（'x='）;
%disp（x）;
y=x;
%%%%%%%%%%%%高斯列主元消元法求解线性方程组Ax=b%%%%%%%%%%%%%%
%A为输入矩阵系数,b为方程组右端系数
%方程组的解保存在x变量中
functiony=Gauss_line（A,b）
formatlong;%设置为长格式显示,显示15位小数
[m,n]=size（A）;
%先检查系数正确性
ifm~=n
error（'矩阵A的行数和列数必须相同'）;
return;
end
ifm~=size（b）
error（'b的大小必须和A的行数或A的列数相同'）;
return;
end
%再检查方程是否存在唯一解
ifrank（A）~=rank（[A,b]）
error（'A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解'）;
return;
end
c=n+1;
A（:
c）=b;%（增广）
fork=1:
n-1
[r,m]=max（abs（A（k:
n,k）））;%选主元
m=m+k-1;%修正操作行的值
if（A（m,k）~=0）
if（m~=k）
A（[km],:
）=A（[mk],:
）;%换行
end
A（k+1:
n,k:
c）=A（k+1:
n,k:
c）-（A（k+1:
n,k）/A（k,k））*A（k,k:
c）;%消去
end
end
x=zeros（length（b）,1）;%回代求解
x（n）=A（n,c）/A（n,n）;
fork=n-1:
-1:
1
x（k）=（A（k,c）-A（k,k+1:
n）*x（k+1:
n））/A（k,k）;
end
y=x;
formatshort;%设置为默认格式显示,显示5位
（2）建立MATLAB界面
利用MATLAB的GUI建立如下界面求解线性方程组：
详见程序。
五、计算实例、数据、结果、分析
下面我们对以上的结果进行测试，求解：
输入数据后点击和，得到如下结果：
更改以上数据进行测试，求解如下方程组：
得到如下结果：
六、实验中遇到的问题及解决办法
在本实验中，遇到的问题主要有两个：
（1）如何将上述的Gauss消元法的算法在MATLAB中实现
针对此问题我借鉴了网上以及课本上的算法的MATLAB实现的程序；
（2）如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解
针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完成了此项任务。
七、实验结论
通过以上的测试，我们发现以上算法和程序能够求出线性方程组的比较精确解。
八、参考文献
[1]杨大地，王开荣.2006.数值分析.北京:
科学出版社
[2]何光辉.2008.数值分析实验.重庆大学数理学院数学实验教学中心
[3]百度文库，百度知道
教师签名
年月日
课程
名称
数值分析
实验项目
名称
插值方法
实验项目类型
验证
演示
综合
设计
其他
指导
教师
何光辉
成绩
√
一、实验目的：
（1）学会拉格朗日插值、牛顿插值等基本方法
（2）设计出相应的算法，编制相应的函数子程序
（3）会用这些函数解决实际问题
二、实验内容
（1）设计拉格朗日插值算法，编制并调试相应的函数子程序
（2）设计牛顿插值算法，编制并调试相应的函数子程序
（3）给定函数四个点的数据如下：
X
1.1
2.3
3.9
5.1
Y
3.887
4.276
4.651
2.117
试用拉格朗日插值确定函数在x=2.101，4.234处的函数值。
（4）已知用牛顿插值公式求的近似值。
三、实验原理
（1）拉格朗日插值
n次拉格朗日插值多项式为：
Ln（x）=y0l0（x）+y1l1（x）+y2l2（x）+…+ynln（x）
n=1时，称为线性插值，L1（x）=y0（x-x1）/（x0-x1）+y1（x-x0）/（x1-x0）=y0+（y1-x0）（x-x0）/（x1-x0）
n=2时，称为二次插值或抛物线插值，精度相对高些
L2（x）=y0（x-x1）（x-x2）/（x0-x1）/（x0-x2）+y1（x-x0）（x-x2）/（x1-x0）/（x1-x2）+y2（x-x0）（x-x1）/（x2-x0）/（x2-x1）
对节点xi（i=0,1,…,n）中任一点xk（0<=k<=n）作一n次多项式lk（xk），使它在该点上取值为1，而在其余点xi（i=0,1,…,k-1,k+1,…,n）上为0，则插值多项式为Ln（x）=y0l0（x）+y1l1（x）+y2l2（x）+…+ynln（x）
上式表明：
n个点xi（i=0,1,…,k-1,k+1,…,n）都是lk（x）的零点。
（2）牛顿插值
插商公式
Newton插值多项式为
四、MATLAB软件实现
（1）分别写出lagrange插值法和Newton插值法的求解函数
%%%%%%%%%%%%%%lagrange插值法求解函数%%%%%%%%%%%%%%%%
%x，y为初始数据，z为插值点
functionz=lagrange（x,y,a）
formatlong;%显示15位
n=length（x）;%取长度
%初始计算
s=0;
%进入公式计算
forj=0:
（n-1）
t=1;
fori=0:
（n-1）
ifi~=j
t=t*（a-x（i+1））/（x（j+1）-x（i+1））;
end
end
s=s+t*y（j+1）;
end
z=s;%显示输出结果
formatshort;
%%%%%%%%%%%%%%%Newton插值法求解函数%%%%%%%%%%%%%%
%x，y为初始数据，z为插值点
functionj=Newton（x,y,z）
n=max（size（x））;
l=1;
a=y
（1）;
B=a;
s=1;%一次因子的乘积,预设为1
dx=y;%差商
fori=1:
n-1
dx0=dx;
forj=1:
n-i
dx（j）=（dx0（j+1）-dx0（j））/（x（i+j）-x（j））;
end
df=dx
（1）;
s=s*（z-x（i））;%一次因子乘积
a=a+s*df;%计算各次Newton插值的值
l=l+1;
B=a;%结果保存在变量B中
end
j=B;
（2）建立界面
利用MATLAB中的GUI编程建立如下界面：
详见程序。
五、计算实例、数据、结果、分析
下面我们对以上的问题进行测试：
输入数据：
计算结果如下：
当x=2.101时，
x=4.234时，
同理可以测试（4）中的的值。
六、实验中遇到的问题及解决办法
在本实验中，遇到的问题主要有两个：
（3）如何将上述的插值的算法在MATLAB中实现
针对此问题我借鉴了网上以及课本上的算法的MATLAB实现的程序；
（4）如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解
针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完成了此项任务。
七、实验结论
通过以上的测试，我们发现以上算法和程序能够求出插值的比较精确解。
八、参考文献
[1]杨大地，王开荣.2006.数值分析.北京:
科学出版社
[2]何光辉.2008.数值分析实验.数理学院数学实验教学中心
[3]百度文库，百度知道
教师签名
年月日
课程
名称
数值分析
实验项目
名称
数值微积分
实验项目类型
验证
演示
综合
设计
其他
指导
教师
何光辉
成绩
√
一、实验目的：
（1）学会复化梯形、复化辛浦生求积公式的应用
（2）设计出相应的算法，编制相应的函数子程序
（3）会用这些函数解决实际问题
二、实验内容
（1）设计复化梯形公式求积算法，编制并调试相应的函数子程序
（2）设计复化辛浦生求积算法，编制并调试相应的函数子程序
（4）分别用复化梯形公式和复化辛浦生公式计算定积分
三、实验原理
（1）复化梯形求积公式
开始
定义；给出
结束
输入
输出
图1复化梯形求积公式算法的流程图
Step1给出被积函数、区间端点和等分数；
Step2求出；
Step3计算；
Step4得
（2）复化辛普森求积公式
开始
定义；给出
结束
输入
输出
图2复化辛普森求积公式算法的流程图
Step1给出被积函数、区间端点和等分数；
Step2求出；
Step3计算；
Step4得
四、MATLAB软件实现
（1）分别写出复化梯形和复化辛浦生求积的求解函数
%%%%%%%%%%%%%%%%%%复化梯形公式求积分值%%%%%%%%%%%%%%%%%
functionT=trap（f,a,b）
%f为积分函数
%[a,b]为积分区间
%n是等分区间份数
n=200;
h=（b-a）/n;%步长
T=0;
fork=1:
（n-1）
x0=a+h*k;
T=T+limit（f,x0）;
end
T=h*（limit（f,a）+limit（f,b））/2+h*T;
T=double（T）;
%%%%%%%%%%%%%Simpson公式求积分值%%%%%%%%%%%%%%%%%%%%%
functionS=simpson（f,a,b）
%f为积分函数
%[a,b]为积分区间
%n是等分区间份数
n=200;
h=（b-a）/（2*n）;%步长
s1=0;
s2=0;
fork=1:
n
x0=a+h*（2*k-1）;
s1=s1+limit（f,x0）;
end
fork=1:
（n-1）
x0=a+h*2*k;
s2=s2+limit（f,x0）;
end
S=h*（limit（f,a）+limit（f,b）+4*s1+2*s2）/3;
S=double（S）;
（2）建立界面
利用MATLAB中的GUI编程建立如下界面：
详见程序。
五、计算实例、数据、结果、分析
下面我们对以上的问题进行测试：
输入数据：
计算结果如下：
六、实验中遇到的问题及解决办法
在本实验中，遇到的问题主要有两个：
（5）如何将上述的积分算法在MATLAB中实现
针对此问题我借鉴了网上以及课本上的算法的MATLAB实现的程序；
（6）如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解
针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完成了此项任务。
七、实验结论
通过以上的测试，我们发现以上算法和程序能够求出积分的比较精确解。
八、参考文献
[1]杨大地，王开荣.2006.数值分析.北京:
科学出版社
[2]何光辉.2008.数值分析实验.数理学院数学实验教学中心
[3]百度文库，百度知道
教师签名
年月日
课程
名称
数值分析
实验项目
名称
常微分方程的数值解法
实验项目类型
验证
演示
综合
设计
其他
指导
教师
何光辉
成绩
√
一、实验目的：
（1）学会欧拉方法和四阶龙格-库塔方法的使用
（2）设计出相应的算法，编制相应的函数子程序
（3）会用这些函数解决实际问题
二、实验内容
用欧拉方法和四阶龙格-库塔方法求解微分方程初值问题：
y’=sin（y）*x，y（0）=10，求y
（1）。
三、实验原理
（1）欧拉法
欧拉法是解初值问题的最简单的数值方法。
从（9.2）式由于y（x0）=y0已给定，因而可以算出
设x1=h充分小，则近似地有：
记
从而我们可以取
作为y（x1）的近似值。
利用y1及f（x1,y1）又可以算出y（x2）的近似值：
一般地，在任意点xn+1=（n+1）h处y（x）的近似值由下式给出
这就是欧拉法的计算公式，h称为步长。
在实际计算时，可将欧拉法与梯形法则相结合，计算公式为：
（2）四阶龙格-库塔方法
四阶龙格-库塔法求解公式如下：
四、MATLAB软件实现
（1）分别写出欧拉方法和四阶龙格-库塔方法求解微分方程的求解函数
%%%%%%%%%Euler法，初值y（a）=c%%%%%%%%%%%%%%%
functiony=Euler（f,a,b,c）
n=1000;
h=（b-a）/n;
X=a:
h:
b;
Y=zeros（1,n+1）;
Y
（1）=c;
fori=2:
n+1
x=X（i-1）;
y=Y（i-1）;
Y（i）=Y（i-1）+eval（f）*h;
end
%%%%%%%%%%%%%%%%%%%ËĽ×Áú¸ñ-¿âËþ·½·¨%%%%%%%%%%%%%%%%
functiony=RK（f,a,b,c）
n=1000;
h=（b-a）/n;
X=a:
h:
b;
Y=zeros（1,n+1）;
Y
（1）=c;
fori=1:
n
x=X（i）;
y=Y（i）;
K1=h*eval（f）;
x=x+h/2;
y=y+K1/2;
K2=h*eval（f）;
x=x;
y=Y（i）+K2/2;
K3=h*eval（f）;
x=X（i）+h;
y=Y（i）+K3;
K4=h*eval（f）;
Y（i+1）=Y（i）+（K1+2*K2+2*K3+K4）/6;
end
（2）建立界面
利用MATLAB中的GUI编程建立如下界面：
详见程序。
五、计算实例、数据、结果、分析
下面我们对以上的问题进行测试：
输入数据：
计算结果如下：
六、实验中遇到的问题及解决办法
在本实验中，遇到的问题主要有两个：
（7）如何将上述的求解微分方程的算法在MATLAB中实现
针对此问题我借鉴了网上以及课本上的算法的MATLAB实现的程序；
（8）如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解
针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完成了此项任务。
七、实验结论
通过以上的测试，我们发现以上算法和程序能够求出微分方程组的比较精确解。
八、参考文献
[1]杨大地，王开荣.2006.数值分析.北京:
科学出版社
[2]何光辉.2008.数值分析实验.数理学院数学实验教学中心
[3]百度文库，百度知道
教师签名
年月日
课程
名称
数值分析
实验项目
名称
估计水塔的水流量
实验项目类型
验证
演示
综合
设计
其他
指导
教师
何光辉
成绩
√
一、实验目的：
（1）学会对实际问题的分析方法
（2）学会利用所学的知识解决实际问题
（3）设计出相应的算法，编制相应的应用程序
二、实验内容
某居民区，其自来水是有一个圆柱形水塔提供，水塔高12.2m，塔的直径为17.4m，水塔是由水泵根据水塔中的水位自动加水，一般水泵每天工作两次。
按照设计，当水塔中的水位降低至最低水位，约8.2m时，水泵自动启动加水。
当水位升至最高水位，约10.8m时，水泵停止工作。
表略。
三、实验原理
计算中将流量定义为单位时间流出的水的高度乘以水塔横截面积。
把时间分成5段：
第1未供水段、水泵开启第1段、第2未供水段、水泵开启第2段、第3未供水段。
先直接对第1、2、3未供水段进行5次曲线拟合。
再对得到的曲线分别求导，取得流速（即单位时间内流出的水的高度）。
水泵开启第1、2段，分别在两端各取两个点，用时刻流速进行拟合得到这两段的流速。
流速乘以水塔横截面积就得到任何时刻的水流量。
对其进行分段积分，求和得到一天的总水流量。
四、MATLAB软件实现
（1）程序
functionpushbutton1_Callback（hObject,eventdata,handles）
%hObjecthandletopushbutton1（seeGCBO）
%eventdatareserved-tobedefinedinafutureversionofMATLAB
%handlesstructurewithhandlesanduserdata（seeGUIDATA）
figure
（1）;
x=[0,3316,6635,10619,13937,17921,21240,25223,28543,32284,39435,43318,46636,49953,53936,57254,60574,64554,68535,71854,75021,85968,89953,93270];
y=[31.75,31.10,30.54,29.94,29.55,28.92,28.50,27.87,27.52,26.97,35.50,34.45,33.50,32.67,31.56,30.81,30.12,29.27,28.42,27.67,26.97,34.75,33.89,33.40];
t=x/3600;%时间单位为小时
h=y/3.281;%水位高度为米
x1=t（1:
10）;
y1=h（1:
10）;
f1=polyfit（x1,y1,5）;
t1=0:
0.01:
t（10）;
h1=polyval（f1,t1）;
plot（x1,y1,'o',t1,h1,'k'）;
xlabel（'时间（h）'）;
ylabel（'水位（m）'）;
title（'第一阶段供水的时间水位图'）
functionpushbutton2_Callback（hObject,eventdata,handles）
%hObjecthandletopushbutton2（seeGCBO）
%eventdatareserved-tobedefinedinafutureversionofMATLAB
%handlesstructurewithhandlesanduserdata（seeGUIDATA）
figure
（2）;
x=[0,3316,6635,10619,13937,17921,21240,25223,28543,32284,39435,43318,46
升级会员 