完整word版计算方法A上机实验报告.docx
《完整word版计算方法A上机实验报告.docx》由会员分享，可在线阅读，更多相关《完整word版计算方法A上机实验报告.docx（27页珍藏版）》请在冰豆网上搜索。
完整word版计算方法A上机实验报告
计算方法A上机实验报告
姓名：
苏福班级：
硕4020学号：
3114161019
一、上机练习目的
1）复习和巩固数值计算方法的基本数学模型，全面掌握运用计算
机进行数值计算的具体过程及相关问题。
2）利用计算机语言独立编写、调试数值计算方法程序，培养学生
利用计算机和所学理论知识分析解决实际问题的能力。
二、上机练习任务
1）利用计算机语言编写并调试一系列数值方法计算通用程序，并
能正确计算给定题目，掌握调试技能。
2）掌握文件使用编程技能，如文件的各类操作，数据格式设计、
通用程序运行过程中文件输入输出运行方式设计等。
3）写出上机练习报告。
三、上机题目
1.共轭梯度法求解线性方程组。
（第三章）
2.三次样条插值（第四章）
3.龙贝格积分（第六章）
4.四阶龙格－库塔法求解常微分方程的初值问题
四、上机报告
题目1：
共轭梯度法求解线性方程组
1．算法原理
共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小值的
问题。
从任意给定的初始点出发，沿一组关于矩阵
共轭的方向进行线性搜索，在无舍入误差的假定下，最多迭代
次（其中
为矩阵
的阶数），就可求得二次函数的极小值，也就求得了线性方程组
的解。
定理：
设
是
阶对称正定矩阵，则
是方程组
的解得充分必要条件是
是二次函数
的极小点，即
共轭梯度法的计算公式：
2.程序框图
3.MATLAB编程实现
（1）编写共轭梯度法求解对称正定矩阵的线性方程组见附录（myge.m）：
functionx=myge（A,b）
输入对称正定矩阵及对应的列向量，初始向量设为0，精度取为
。
函数的输出即为由共轭梯度发求解的近似解。
（2）编写具体算例求解（example.m）：
clc
clearall
%例题3.4.2
A0=[2,0,1;0,1,0;1,0,2];
b0=[3,1,3]';
myge（A0,b0）;
%习题3.2
n=100;%矩阵阶数
A=zeros（n,n）;
b=zeros（n,1）;
b
（1）=-1;b（n）=-1;
A（1,1）=-2;A（1,2）=1;A（n,n-1）=1;A（n,n）=-2;
fori=2:
n-1
A（i,i-1）=1;
A（i,i）=-2;
A（i,i+1）=1;
end
myge（A,b）;
算例1（课本例题3.4.2）：
算例2（课后习题计算实习3.2）：
4.算例结果
算例1：
x=
1.000000000000000
1.000000000000000
1.000000000000000
迭代次数：
k=
2
算例2（n=100）：
x=
0.999999999999999
1.000000000000005
0.999999999999982
1.000000000000038
0.999999999999976
1.000000000000013
1.000000000000000
1.000000000000008
0.999999999999984
1.000000000000017
0.999999999999977
1.000000000000052
0.999999999999942
1.000000000000058
0.999999999999974
1.000000000000011
0.999999999999992
1.000000000000028
0.999999999999983
1.000000000000008
1.000000000000021
0.999999999999977
1.000000000000027
0.999999999999987
1.000000000000031
0.999999999999966
1.000000000000043
0.999999999999985
1.000000000000019
0.999999999999993
1.000000000000022
0.999999999999990
1.000000000000019
1.000000000000006
1.000000000000004
1.000000000000011
0.999999999999996
1.000000000000022
0.999999999999991
1.000000000000026
0.999999999999996
1.000000000000018
1.000000000000004
1.000000000000014
1.000000000000012
1.000000000000008
1.000000000000009
1.000000000000015
1.000000000000010
1.000000000000012
1.000000000000012
1.000000000000010
1.000000000000015
1.000000000000007
1.000000000000015
1.000000000000007
1.000000000000012
1.000000000000007
1.000000000000019
1.000000000000001
1.000000000000023
0.999999999999992
1.000000000000021
0.999999999999997
1.000000000000021
0.999999999999984
1.000000000000029
0.999999999999985
1.000000000000030
0.999999999999986
1.000000000000026
0.999999999999974
1.000000000000054
0.999999999999951
1.000000000000062
0.999999999999943
1.000000000000078
0.999999999999934
1.000000000000060
0.999999999999968
1.000000000000021
0.999999999999990
1.000000000000019
0.999999999999999
1.000000000000012
0.999999999999978
1.000000000000042
0.999999999999963
1.000000000000057
0.999999999999933
1.000000000000063
0.999999999999960
1.000000000000034
0.999999999999974
1.000000000000025
0.999999999999969
1.000000000000054
0.999999999999955
1.000000000000022
0.999999999999988
迭代次数：
k=
50
对于n=200,400，编写的程序的计算结果显示可靠，满足精度要求。
题目2:
三次样条插值
1.算法原理
设在区间
上给定
个节点
，在节点
处的函数值
为
。
若函数
满足如下三条：
（1）在每个子区间
上，
是三次多项式；
（2）
；
（3）在区间
上，
的二阶导数
连续。
则称
为函数
在区间
上的三次样条插值函数。
子区间
上的
的表达式为：
关于参数
的方程组（三弯矩方程组）：
牛顿插值多项式：
构造牛顿插值多项式首先列出差商表，进而由差商表写出牛顿插值多项式。
阶差商为：
零阶差商和一阶差商：
2.程序框图
3.MATLAB编程实现
给定函数
，取等距节点。
编写三次样条插值函数见附录（myspline.m）：
functionSx=myspline（n）
输入为区间个数，输出为各区间上的三次样条插值函数。
编写牛顿插值函数见附录（myNewton.m）：
functionNx=myNewton（n）
输入为区间个数，输出为牛顿插值函数。
编写main函数（main.m）：
clc
clearall
x=-1:
0.001:
1;
y=1./（1+25*x.^2）;
plot（x,y,'k-'）;holdon;
Nx5=subs（myNewton（5）,x）;
plot（x,Nx5,'c-.'）;holdon;
Nx10=subs（myNewton（10）,x）;
plot（x,Nx10,'r--'）;holdon;
myspline（10）;holdon
legend（'f（x）=1/（1+25x^2）','f（x）=N5（x）','f（x）=N10（x）','f（x））=S10（x）'）
4.算例结果
n=
5
次牛顿插值多项式——
Nx=
0.1538*x+1.058*（x+1.0）*（x+0.6）-1.923*（x+1.0）*（x+0.6）*（x+0.2）+1.202*（x+1.0）*（x+0.6）*（x+0.2）*（x-0.2）+0.1923
n=
10
次牛顿插值多项式——
Nx=
0.1018*x+0.2602*（x+1.0）*（x+0.8）+0.7919*（x+1.0）*（x+0.6）*（x+0.8）+2.687*（x+1.0）*（x+0.6）*（x+0.4）*（x+0.8）-6.363*（x+1.0）*（x+0.6）*（x+0.2）*（x+0.4）*（x+0.8）-17.68*x*（x+1.0）*（x+0.6）*（x+0.2）*（x+0.4）*（x+0.8）+84.84*x*（x+1.0）*（x+0.6）*（x+0.2）*（x+0.4）*（x+0.8）*（x-0.2）-167.9*x*（x+1.0）*（x+0.6）*（x+0.2）*（x+0.4）*（x+0.8）*（x-0.2）*（x-0.4）+220.9*x*（x+1.0）*（x+0.6）*（x+0.2）*（x+0.4）*（x+0.8）*（x-0.2）*（x-0.4）*（x-0.6）-220.9*x*（x+1.0）*（x+0.6）*（x+0.2）*（x+0.4）*（x+0.8）*（x-0.2）*（x-0.4）*（x-0.8）*（x-0.6）+0.1403
n=
10
次三次样条插值多项式——
Sx=
0.08814*x+0.3417*（x+1.0）^3+0.1266
0.1702*x-0.3417*（x+0.6）^3+1.235*（x+0.8）^3+0.1922
0.4665*x+2.071*（x+0.6）^3-1.235*（x+0.4）^3+0.37
0.9637*x-2.071*（x+0.2）^3+15.48*（x+0.4）^3+0.5689
4.679*x-38.99*（x+0.2）^3-15.48*x^3+1.312
38.99*（x-0.2）^3-4.679*x+15.48*x^3+1.312
2.071*（x-0.2）^3-0.9637*x-15.48*（x-0.4）^3+0.5689
1.235*（x-0.4）^3-0.4665*x-2.071*（x-0.6）^3+0.37
0.3417*（x-0.6）^3-1.235*（x-0.8）^3-0.1702*x+0.1922
0.1266-0.3417*（x-1.0）^3-0.08814*x
三种插值多项式与
的曲线在同一坐标系:
题目3：
龙贝格积分
1．算法原理
通过逐步缩小步长
的方法，以使积分的近似值满足精度要求。
在实际中，通过将区间
间逐次对半分来实现。
在变步长积分的思想上，对于复化梯形积分：
对于复化辛普森求积公式：
同理可以证明：
再令：
即为龙贝格积分公式。
龙贝格积分是将区间
逐次分半，逐次递推计算，以得到较高精度的积分近似值。
龙贝格的积分公式如下：
若
2.程序框图
3.MATLAB编程实现
（1）编写自定义函数（f.m）：
functionfun=f（x）
fun=log（1+x）/x;%改变函数表达式求不同的函数的积分
取习题6.2的积分函数为例采用龙贝格积分方法求积分。
（2）编写龙贝格积分程序见附录（myromberg.m）：
4.算例结果
（1）
R_k=
0.693147978563680
（2）
R_k=
0.272198867426936
（3）
R_k=
0.822467764886147
（4）
R_k=
1.370762774269822
题目4：
四阶龙格—库塔法求解常微分方程的初值问题
1.算法原理
对于高阶常微分方程的初值问题：
将其转化为一阶微分方程组求解：
标准的四级四阶龙格—库塔法的向量形式：
其分量形式：
2.程序框图
3.MATLAB编程实现
编写4级4阶龙格—库塔积分法求解（步长0.05）：
转化为一阶常微分方程组：
编写MATLAB程序myRK.m（见附录），求出总长度为5的数值积分结果，并与解析解
作比较。
4.算例结果
龙格-库塔积分结果与实际函数值比较。
五、上机总结
通过本次计算方法A课程的上机实验，复习和巩固数值计算方法的基本数学模型，对共轭梯度法求解线性方程组，三次样条插值，牛顿插值，龙贝格积分法，标准四级四阶龙格——库塔法求解高阶常微分方程的初值问题等知识有深刻的理解与认识，并通过MATLAB编程实现各种数值算法，对数值分析的应用有清晰的认识。
附录
1.共轭梯度法myge.m
%共轭梯度法
functionx=myge（A,b）
length=size（b,1）;
x0=zeros（length,1）;%迭代初始向量
r0=b-A*x0;
d0=r0;
fork=1:
length
a=r0'*r0/（d0'*A*d0）;
x0=x0+a*d0;
r1=b-A*x0;
ifnorm（r1）<10^（-8）
x=x0;
formatlong
x
disp（'迭代次数'）;
k
break;
end
p=norm（r1）^2/norm（r0）^2;
d0=r1+p*d0;
r0=r1;
end
2.三次样条插值myspline.m
%三次样条插值
functionSx=myspline（n）
a=-1;b=1;
xi=a:
（b-a）/n:
b;
yi=1./（1+25*xi.^2）;
fori=1:
n
h（i）=xi（i+1）-xi（i）;
end
fori=1:
（n-1）
u（i）=h（i）/（h（i）+h（i+1））;
v（i）=1-u（i）;
end
fori=2:
n
d（i-1）=6*（（yi（i+1）-yi（i））/h（i）-（yi（i）-yi（i-1））/h（i-1））/（h（i）+h（i-1））;
end
zu=zeros（n-1,1）;zy=zeros（n-1,1）;M=zeros（n+1,1）;
M（n+1,1）=0;M（1,1）=0;%第一种边界条件（设为自然边界条件）
d
（1）=d
（1）-u
（1）*M（1,1）;d（n-1）=d（n-1）-v（n-1）*M（n+1,1）;
zu（1,1）=2;zy（1,1）=d
（1）;
%追赶法求解三弯矩方程组
fori=2:
（n-1）
l=u（i）/zu（i-1,1）;
zu（i,1）=2-l*v（i-1）;
zy（i,1）=d（i）-l*zy（i-1,1）;
l=0;
end
M（n,1）=zy（n-1,1）/zu（n-1,1）;
forj=1:
（n-2）
M（n-j,1）=（zy（n-j-1,1）-v（n-j-1）*M（n-j+1,1））/zu（n-j-1,1）;
end
M;
symsx
Sxi=0;
set_Sx=[];
forj=1:
n
Sxi=（xi（j+1）-x）^3/（6*h（j））*M（j,1）+（x-xi（j））^3/（6*h（j））*M（j+1,1）+（yi（j）-h（j）^2/6*M（j,1））*（xi（j+1）-x）/h（j）+（yi（j+1）-h（j）^2/6*M（j+1,1））*（x-xi（j））/h（j）;
set_Sx=[set_Sx;Sxi];
end
set_Sx;
Sx=set_Sx;
forj=1:
n
x_Sx=xi（j）:
0.0001:
xi（j+1）;
y_Sx=subs（set_Sx（j）,x_Sx）;
plot（x_Sx,y_Sx,'b:
'）;
holdon
end
n
disp（'次三次样条插值多项式——'）
Sx
Vpa（Sx,4）
3.牛顿插值myNewton.m
%myNewton.m
functionNx=myNewton（n）
a=-1;b=1;
xi=a:
（b-a）/n:
b;
yi=1./（1+25*xi.^2）;
table=zeros（n+1,n+1）;
table（:
1）=yi';%差商表table
forj=2:
n+1
fori=j:
n+1
table（i,j）=（table（i,j-1）-table（i-1,j-1））/（xi（i）-xi（i-j+1））;
end
end
symsx
Nx=table（1,1）;
Tx=1;
fork=2:
n+1
Tx=Tx*（x-xi（k-1））;
Nx=Nx+table（k,k）*Tx;
end
n
disp（'次牛顿插值多项式——'）
Nx
vpa（Nx,4）
4.龙贝格积分myromberg.m
%myromberg.m
clc
clearall
%积分区间
a=0;b=1;
a=a+10^（-10）;%a取极小值以使得f（0）=1，避免sin（0）/o等
f（a）;
T_1=（b-a）*（f（a）+f（b））/2;
T_2=T_1/2+（b-a）/2*f（a+（b-a）/2）;
T_4=T_2/2+（b-a）/4*（（f（a+（b-a）/4））+f（a+（b-a）*3/4））;
T_8=T_4/2+（b-a）/8*（（f（a+（b-a）/8））+f（a+（b-a）*3/8）+f（a+（b-a）*5/8）+a+（b-a）*7/8）;
S_1=T_2+1/（4-1）*（T_2-T_1）;S_2=T_4+1/（4-1）*（T_4-T_2）;S_4=T_8+1/（4-1）*（T_8-T_4）;
C_1=S_2+1/（4^2-1）*（S_2-S_1）;C_2=S_4+1/（4^2-1）*（S_4-S_2）;
R_1=C_2+1/（4^3-1）*（C_2-C_1）;
Tsum=0;
R_k=R_1+0.2;
k=3;
%开始迭代
whileabs（R_k-R_1）>1e-6%精度设为1e-6
R_1=R_k;
fori=1:
2^k
Tsum=Tsum+f（a+（2*i-1）*（b-a）/2^（k+1））;
end
T_k=T_8/2+（b-a）/2^（k+1）*Tsum;
S_k=T_k+1/（4-1）*（T_k-T_8）;
C_k=S_k+1/（4^2-1）*（S_k-S_4）;
R_k=C_k+1/（4^3-1）*（C_k-C_2）;
T_8=T_k;
S_4=S_k;
C_2=C_k;
Tsum=0;
k=k+1;
end
k
formatlong
R_k
5.四阶龙格—库塔法myRK.m
%myRK.m
clc
clearall
h=0.05;
n=5/h;
y=zeros（3,n）;
y（:
1）=[-1;3;2];
K1=zeros（3,1）;K2=zeros（3,1）;K3=zeros（3,1）;K4=zeros（3,1）;
fori=1:
n
x=i*h;
forj=1:
3
K1
（1）=h*y（2,i）;K1
（2）=h*y（3,i）;K1（3）=h*（y（3,i）+y（2,i）-y（1,i）+2*x-3）;
K2
（1）=h*（y（2,i）+K1
（2）/2）;
K2
（2）=h*（y（3,i）+K1（3）/2）;
K2（3）=h*（y（3,i）+K1（3）/2+y（2,i）+K1
（2）/2-（y（1,i）+K1
（1）/2）+2*（x+h/2）-3）;
K3
（1）=h*（y（2,i）+K2
（2）/2）;
K3
（2）=h*（y（3,i）+K2（3）/2）;
K3（3）=h*（y（3,i）+K2（3）/2+y（2,i）+K2
（2）/2-（y（1,i）+K2
（1）/2）+2*（x+h/2）-3）;
K4
（1）=h*（y（2,i）+K3
（2））;
K4
（2）=h*（y（3,i）+K3（3））;
K4（3）=h*（y（3,i）+K3（3）+y（2,i）+K3
（2）-（y（1,i）+K3
（1））+2*（x+h）-3）;
end
fork=1:
3
y（k,i+1）=y（k,i）+（K1（k）+2*K2（k）+2*K3（k）+K4（k））/6;
end
end
xh=h:
h:
（n+1）*h;
plot（xh,y（1,:
）,'r.'）;holdon
xi=h:
0.01:
（n+1）*h;
yi=xi.*exp（xi）+2*xi-1;
plot（xi,yi,'k'）;
holdon
legend（'4级R-K','y（x）=xe^x+2x-1'）
升级会员 