数值分析上机题目.docx
《数值分析上机题目.docx》由会员分享，可在线阅读，更多相关《数值分析上机题目.docx（24页珍藏版）》请在冰豆网上搜索。
数值分析上机题目
--本页仅作为文档封面，使用时请直接删除即可----内页可以根据需求调整合适字体及大小--
数值分析上机题目4（总21页）
实验一
实验项目：
共轭梯度法求解对称正定的线性方程组
实验内容：
用共轭梯度法求解下面方程组
（1）
迭代20次或满足
时停止计算。
编制程序：
储存m文件
function[x,k]=CGmethod（A,b）
n=length（A）;x=2*ones（n,1）;r=b-A*x;rho=r'*r;
k=0;
whilerho>10^（-11）&k<1000
k=k+1;
ifk==1
p=r;
else
beta=rho/rho1;
p=r+beta*p;
end
w=A*p;
alpha=rho/（p'*w）;
x=x+alpha*p;
r=r-alpha*w;
rho1=rho;
rho=r'*r;
end
运行程序：
clear,clc
A=[2-100;-13-10;0-14-1;00-15];
b=[3-215]';
[x,k]=CGmethod（A,b）
运行结果：
x=
（2）
A是1000阶的Hilbert矩阵或如下的三对角矩阵，
A[i,i]=4，A[i,i-1]=A[i-1,i]=-1，i=2,3,..,n
b[1]=3,b[n]=3,b[i]=2，i=2,3,…,n-1
迭代10000次或满足
时停止计算。
编制程序：
储存m文件
function[x,k]=CGmethod_1（A,b）
n=length（A）;x（1:
n,1）=0;r=b-A*x;r1=r;
k=0;
whilenorm（r1,1）>10^（-7）&k<10^4
k=k+1;
ifk==1
p=r;
else
beta=（r1'*r1）/（r'*r）;p=r1+beta*p;
end
r=r1;
w=A*p;
alpha=（r'*r）/（p'*w）;
x=x+alpha*p;
r1=r-alpha*w;
end
运行程序：
clear,clc
n=1000;
A=hilb（n）;
b=sum（A'）';
[x,k]=CGmethod（A,b）
实验二
1、实验目的：
用复化Simpson方法、自适应复化梯形方法和Romberg方法求数值积分。
实验内容：
计算下列定积分
（1）
（2）
（3）
实验要求：
（1）分别用复化Simpson公式、自适应复化梯形公式计算要求绝对误差限为
，输出每种方法所需的节点数和积分近似值，对于自适应方法，显示实际计算节点上离散函数值的分布图；
（2）分析比较计算结果。
程序：
symsx
f=x^6/10-x^2+x%定义函数f（x）
n=input（'输入所求导数阶数:
'）
f2=diff（f,x,n）%求f（x）的n阶导数
（1）复化梯形
clc
clear
symsx%定义自变量x
f=inline（'x^6/10-x^2+x','x'）%定义函数f（x）=x*exp（x），换函数时只需换该函数表达式即可
f2=inline（'3*x^4-2','x'）%定义f（x）的二阶导数,输入程序1里求出的f2即可。
f3='-（3*x^4-2）'%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，以便求最大值
e=*10^（-7）%精度要求值
a=0%积分下限
b=2%积分上限
x1=fminbnd（f3,1,2）%求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值
forn=2:
1000000%求等分数n
Rn=-（b-a）/12*（（b-a）/n）^2*f2（x1）%计算余项
ifabs（Rn）break%符合要求时结束
end
end
h=（b-a）/n%求h
Tn1=0
fork=1:
n-1%求连加和
xk=a+k*h
Tn1=Tn1+f（xk）
end
Tn=h/2*（（f（a）+2*Tn1+f（b）））
z=exp
（2）
R=Tn-z%求已知值与计算值的差
stem（xk,Tn1）;
fprintf（'用复化梯形算法计算的结果Tn='）
disp（Tn）
fprintf（'等分数n='）
disp（n）%输出等分数
fprintf（'已知值与计算值的误差R='）
disp（R）
（2）复化Simpson
clc
clear
symsx%定义自变量x
f=inline（'x^6/10-x^2+x','x'）%定义函数f（x）=x*exp（x），换函数时只需换该函数表达式即可
f2=inline（'36*x^2','x'）%定义f（x）的四阶导数，输入程序1里求出的f2即可
f3='-（36*x^2）'%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，一边求最大值
e=5*10^（-8）%精度要求值
a=0%积分下限
b=2%积分上限
x1=fminbnd（f3,1,2）%求负的四阶导数的最小值点，也就是求四阶导数的最大值点对应的x值
forn=2:
1000000%求等分数n
Rn=-（b-a）/180*（（b-a）/（2*n））^4*f2（x1）%计算余项
ifabs（Rn）break%符合要求时结束
end
end
h=（b-a）/n%求h
Sn1=0
Sn2=0
fork=0:
n-1%求两组连加和
xk=a+k*h
xk1=xk+h/2
Sn1=Sn1+f（xk1）
Sn2=Sn2+f（xk）
end
Sn=h/6*（f（a）+4*Sn1+2*（Sn2-f（a））+f（b））%因Sn2多加了k=0时的值，故减去f（a）
z=exp
（2）
R=Sn-z%求已知值与计算值的差
fprintf（'用Simpson公式计算的结果Sn='）
disp（Sn）
fprintf（'等分数n='）
disp（n）
fprintf（'已知值与计算值的误差R='）
disp（R）
结果：
用复化梯形算法计算的结果Tn=
等分数n=24764
已知值与计算值的误差R=
用Simpson公式计算的结果Sn=
等分数n=76
已知值与计算值的误差R=
用复化梯形算法计算的结果Tn=
等分数n=1119
已知值与计算值的误差R=
用Simpson公式计算的结果Sn=
等分数n=8
已知值与计算值的误差R=
用复化梯形算法计算的结果Tn=
等分数n=1000000
已知值与计算值的误差R=
用Simpson公式计算的结果Sn=
等分数n=10647
已知值与计算值的误差R=
分析：
在处理问题时，复化Simpson要比复化梯度计算速度要快很多。
2、实验目的：
高斯数值积分方法用于积分方程求解。
实验内容：
线性的积分方程的数值求解，可以被转化为线性代数方程组的求解问题。
而线性代数方程组所含未知数的个数，与用来离散积分的数值方法的节点个数相同。
在节点数相同的前提下，高斯数值积分方法有较高的代数精度，用它通常会得到较好的结果。
对第二类Fredholm积分方程
首先将积分区间[a,b]等分成n份，在每个子区间上离散方程中的积分就得到线性代数方程组。
实验要求：
分别使用如下方法，离散积分方程中的积分
1.复化梯形方法；2.复化辛甫森方法；3.复化高斯方法。
求解如下的积分方程
，方程的准确解为
，
并比较各算法的优劣。
程序结果：
当迭代次数n=1时
精确解
复化梯形方法
复化辛甫森方法
复化高斯方法
复化梯形方法的平均误差err=
复化辛甫森方法的平均误差err=
复化高斯方法的平均误差err=
当迭代次数n=5时，
精确解
复化梯形方法
复化辛甫森方法
复化高斯方法
复化梯形方法的平均误差err=
复化辛甫森方法和复化高斯方法的平均误差err=0
可以看出，复化高斯方法得到的结果精度最高，复化辛普森方法比复化梯形方法的精度要高，
程序：
clear,clc
a=0;b=1;
n=5;
figure
fun1（a,b,n）;
holdon
a=0;b=1;
n=5;
figure
fun2（a,b,n）;
holdon
figure
fun3（a,b,n）;
编制m函数：
functiony=Kfun（t,s）
y=2/（exp
（1）-1）*exp（t）;
functiony=ffun（t）
y=-exp（t）;
functiony=Fexc（t）
%精确解
y=exp（t）;
function[x,w]=fhgauss（a,b,n）
h=（b-a）/n;
x1=a:
h:
b;
x=zeros（1,2*n）;
w=x;
fori=1:
n
[x（2*i-1:
2*i）,w（2*i-1:
2*i）]=GaussLegendre（x1（i）,x1（i+1）,2）;
end
function[x,w]=fhsimpson（a,b,n）
h=（b-a）/n;
x=a:
h/2:
b;
w=x;
w
（1）=h/6;
w（2*n+1）=h/6;
fori=1:
n
w（2*i）=2/3*h;
ifiw（2*i+1）=h/3;
end
end
function[x,w]=fhtrapz（a,b,n）
h=（b-a）/n;
x=a:
h:
b;
fori=1:
n+1
ifi==1||i==n+1
w（i）=h/2;
else
w（i）=h;
end
end
function[x,w]=GaussLegendre（a,b,n）
i=1:
n-1;
c=i./sqrt（4*i.^2-1）;
CM=diag（c,1）+diag（c,-1）;
[VL]=eig（CM）;
[xind]=sort（diag（L））;
V=V（:
ind）';
w=2*V（:
1）.^2;
x=x*（b-a）/2+（b+a）/2;
w=w*（b-a）/2;
functionfun1（a,b,n）
[x1,w1]=fhtrapz（a,b,n）;
n1=4;
n=n+1;
h=（b-a）/n1;
x=a:
h:
b;
y0=Fexc（x）;
A=zeros（n,n）;
B=zeros（n,1）;
fori=1:
n
B（i）=ffun（x1（i））;
end
fori=1:
n
forj=1:
n
A（i,j）=w1（j）*Kfun（x1（i）,x1（j））;
end
end
A=eye（n）-A;
y1=（A\B）';
yN=x;
fori=1:
n1+1
yN（i）=ffun（x（i））;
forj=1:
n
yN（i）=yN（i）+w1（j）*Kfun（x（i）,x1（j））*y1（j）;
end
end
fprintf（'数值解:
\n'）
yN
fprintf（'精确解:
\n'）
y0
plot（x,yN,'x',x,y0,'.'）;
h=legend（'复化值','真实值'）;
functionfun2（a,b,n）
[x1,w1]=fhsimpson（a,b,n）;
n=2*n+1;
n1=50;
h=（b-a）/n1;
x=a:
h:
b;
y0=Fexc（x）;
A=zeros（n,n）;
B=zeros（n,1）;
fori=1:
n
B（i）=ffun（x1（i））;
end
fori=1:
n
forj=1:
n
A（i,j）=w1（j）*Kfun（x1（i）,x1（j））;
end
end
A=eye（n）-A;
y1=（A\B）';
yN=x;
fori=1:
n1+1
yN（i）=ffun（x（i））;
forj=1:
n
yN（i）=yN（i）+w1（j）*Kfun（x（i）,x1（j））*y1（j）;
end
end
fprintf（'数值解:
\n'）
yN
fprintf（'精确解:
\n'）
y0
plot（x,yN,'x',x,y0,'.'）;
h=legend（'复化值','真实值'）;
functionfun3（a,b,n）
[x1,w1]=fhgauss（a,b,n）;
n=2*n;
n1=4;
h=（b-a）/n1;
x=a:
h:
b;
y0=Fexc（x）;
A=zeros（n,n）;
B=zeros（n,1）;
fori=1:
n
B（i）=ffun（x1（i））;
end
fori=1:
n
forj=1:
n
A（i,j）=w1（j）*Kfun（x1（i）,x1（j））;
end
end
A=eye（n）-A;
y1=（A\B）';
yN=x;
fori=1:
n1+1
yN（i）=ffun（x（i））;
forj=1:
n
yN（i）=yN（i）+w1（j）*Kfun（x（i）,x1（j））*y1（j）;
end
end
fprintf（'数值解:
\n'）
yN
fprintf（'精确解:
\n'）
y0
plot（x,yN,'x',x,y0,'.'）;
h=legend（'复化值','真实值'）;
图一
图二
图三
实验三
1、对常微分方程初值问题
分别使用Euler显示方法、改进的Euler方法和经典RK法和四阶Adams预测-校正算法，求解常微分方程数值解，并与其精确解
进行作图比较。
其中多步法需要的初值由经典RK法提供。
程序：
子程序
functionE=Euler（fun,x0,y0,h,N）
x=zeros（1,N+1）;
y=zeros（1,N+1）;
x
（1）=x0;
y
（1）=y0;
forn=1:
N
x（n+1）=x（n）+h;
y（n+1）=y（n）+h*feval（fun,x（n）,y（n））;
end
x1=x
y1=y
plot（x1,y1,'k'）
holdon
functionE=Euler_modify（fun,x0,y0,h,N）
x=zeros（1,N+1）;
y=zeros（1,N+1）;
x
（1）=x0;y
（1）=y0;
forn=1:
N
x（n+1）=x（n）+h;
z0=y（n）+h*feval（fun,x（n）,y（n））;
y（n+1）=y（n）+h/2*（feval（fun,x（n）,y（n））+feval（fun,x（n+1）,z0））;
end
x2=x
y2=y
plot（x2,y2,'g'）
function[x,y]=Rk_N4（f,x0,y0,h,N）
x=zeros（1,N+1）;
y=zeros（1,N+1）;
x
（1）=x0;y
（1）=y0;
forn=1:
N
x（n+1）=x（n）+h;
k1=h*feval（f,x（n）,y（n））;
k2=h*feval（f,x（n）+1/2*h,y（n）+1/2*k1）;
k3=h*feval（f,x（n）+1/2*h,y（n）+1/2*k2）;
k4=h*feval（f,x（n）+h,y（n）+k3）;
y（n+1）=y（n）+1/6*（k1+2*k2+2*k3+k4）;
end
运行以下程序
clc
clear
Euler（'fun',0,1,,100）
Euler_modify（'fun',0,1,,100）
[x,y]=Rk_N4（'fun',0,1,,100）
plot（x,y,'-*r'）
holdon
x1=0:
:
pi;
y1=cos（x1）+sin（x1）
plot（x1,y1,''）
title（'误差分析'）;
xlabel（'x轴'）;
ylabel（'y轴'）;
legend（'Euler','Euler改进法','R_K法','精确'）;
axis（[0pi]）;
gridon
画出图形进行比较：
2、实验目的：
Lorenz问题与混沌
实验内容：
考虑著名的Lorenz方程
其中s,r,b为变化区域有一定限制的实参数。
该方程形式简单，表面上看并无惊人之处，但由该方程揭示出的许多现象，促使“混沌”成为数学研究的崭新领域，在实际应用中也产生了巨大的影响。
实验方法：
先取定初值Y0=（x,y,z）=（0,0,0），参数s=10,r=28,b=8/3，用MATLAB编程数值求解，并与MATLAB函数ods45的计算结果进行对比。
实验要求：
（1）对目前取定的参数值s,r和b，选取不同的初值Y0进行运算，观察计算的结果有什么特点解的曲线是否有界解的曲线是不是周期的或趋于某个固定点
（2）在问题允许的范围内适当改变其中的参数值s,r,b，再选取不同的初始值Y0进行试算，观察并记录计算的结果有什么特点是否发现什么不同的现象
程序：
blzeq函数程序
functionydot=blzeq（t,y）
globalSIGMARHOBETA
A=[-BETA0y
（2）;0-SIGMASIGMA;-y
（2）RHO-1];
ydot=A*y;
主程序：
clf
clc
globalSIGMARHOBETA
SIGMA=10.;
RHO=28.;
BETA=8/3;
axis（[050-3030-3030]）
view（3）
holdon
title（'LorenzAttractor'）
y0=[50,50,50];tfinal=100;
[t,y]=ode23（'blzeq',[0tfinal],y0）;
plot3（y（:
1）,y（:
2）,y（:
3））
（1）实验结果及其分析：
题中的方程与程序中的方程的关系是变量进行了轮换，x换成了y，y换成了z，z换成了x。
原点为原方程的一个奇点，当初始位置稍稍偏离原点如取为[0，eps，0，]（按原方程中的顺序，下同）得到的图像如下：
分析：
这是一个典型的奇怪吸引子的图像，曲线有界，但他不收敛于某一点也不是周期的，而是在两个位置附近来回的跳跃。
（2）取初始位置分别为[20，20，20]，[60，60，60]得到的图像如下：
[20，20，20][60，60，60]
分析：
初始变量值相同时，曲线总是被吸引回奇怪吸引子附近作来回跳跃。
只变化b
取s=10,r=20,b=10初值：
[20，20，20]的图像：
取s=10,r=20,b=5初值：
[20，20，20]的图像：
取s=10,r=20,b=8/9初值：
[20，20，20]的图像：
改变s,r也会有发生类似的情况。
升级会员 