计算方法A上机实验报告.doc

计算方法A上机实验报告
姓名：苏福 班级：硕4020 学号：3114161019

一、上机练习目的
1）复习和巩固数值计算方法的基本数学模型，全面掌握运用计算
机进行数值计算的具体过程及相关问题。
2）利用计算机语言独立编写、调试数值计算方法程序，培养学生
利用计算机和所学理论知识分析解决实际问题的能力。
二、上机练习任务
1）利用计算机语言编写并调试一系列数值方法计算通用程序，并
能正确计算给定题目，掌握调试技能。
2） 掌握文件使用编程技能，如文件的各类操作，数据格式设计、
通用程序运行过程中文件输入输出运行方式设计等。
3）写出上机练习报告。
三、上机题目
1. 共轭梯度法求解线性方程组。 （第三章）
2. 三次样条插值（第四章）
3. 龙贝格积分（第六章）
4. 四阶龙格－库塔法求解常微分方程的初值问题
四、上机报告
题目1：共轭梯度法求解线性方程组
1． 算法原理
共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小值的
问题。从任意给定的初始点出发，沿一组关于矩阵共轭的方向进行线性搜索，在无舍入误差的假定下，最多迭代次（其中为矩阵的阶数），就可求得二次函数的极小值，也就求得了线性方程组的解。
定理：设是阶对称正定矩阵，则是方程组的解得充分必要条件是是二次函数的极小点，即

共轭梯度法的计算公式：

2. 程序框图

3. MATLAB编程实现
（1）编写共轭梯度法求解对称正定矩阵的线性方程组见附录（myge.m）：
function x=myge(A,b)
输入对称正定矩阵及对应的列向量，初始向量设为0，精度取为。函数的输出即为由共轭梯度发求解的近似解。
（2）编写具体算例求解（example.m）：
clc
clear all
%例题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;
for i=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）：
function Sx=myspline(n)
输入为区间个数，输出为各区间上的三次样条插值函数。
编写牛顿插值函数见附录(myNewton.m)：
function Nx=myNewton(n)
输入为区间个数，输出为牛顿插值函数。
编写main函数(main.m)：
clc
clear all
x=-1:0.001:1;
y=1./(1+25*x.^2);
plot(x,y,'k-');hold on;
Nx5=subs(myNewton(5),x);
plot(x,Nx5,'c-.');hold on;
Nx10=subs(myNewton(10),x);
plot(x,Nx10,'r--');hold on;
myspline(10);hold on
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）：
function fun=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
%共轭梯度法
function x=myge(A,b)
length=size(b,1);
x0=zeros(length,1); %迭代初始向量
r0=b-A*x0;
d0=r0;
for k=1:length
a=r0'*r0/(d0'*A*d0);
x0=x0+a*d0;
r1=b-A*x0;
if norm(r1)<10^(-8)
x=x0;
format long
x
disp('迭代次数');
k
break;
end
p=norm(r1)^2/norm(r0)^2;
d0=r1+p*d0;
r0=r1;
end

2. 三次样条插值myspline.m
%三次样条插值
function Sx=myspline(n)
a=-1;b=1;
xi=a:(b-a)/n:b;
yi=1./(1+25*xi.^2);
for i=1:n
h(i)=xi(i+1)-xi(i);
end
for i=1:(n-1)
u(i)=h(i)/(h(i)+h(i+1));
v(i)=1-u(i);
end
for i=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);
%追赶法求解三弯矩方程组
for i=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);
for j=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;
syms x
Sxi=0;
set_Sx=[];
for j=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;
for j=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:');
hold on
end
n
disp('次三次样条插值多项式——')
Sx
Vpa(Sx,4)

3. 牛顿插值myNewton.m
%myNewton.m
function Nx=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
for j=2:n+1
for i=j:n+1
table(i,j)=(table(i,j-1)-table(i-1,j-1))/(xi(i)-xi(i-j+1));
end
end
syms x
Nx=table(1,1);
Tx=1;
for k=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
clear all
%积分区间
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;
%开始迭代
while abs(R_k-R_1)>1e-6 %精度设为1e-6
R_1=R_k;
for i=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
format long
R_k

5. 四阶龙格—库塔法myRK.m
%myRK.m
clc
clear all
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);
for i=1:n
x=i*h;
for j=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
for k=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.');hold on
xi=h:0.01:(n+1)*h;
yi=xi.*exp(xi)+2*xi-1;
plot(xi,yi,'k');
hold on
legend('4级R-K','y(x)=xe^x+2x-1')