数值分析上机第四次作业.docx
《数值分析上机第四次作业.docx》由会员分享，可在线阅读，更多相关《数值分析上机第四次作业.docx（15页珍藏版）》请在冰豆网上搜索。
数值分析上机第四次作业
数值分析上机第四次作业
数值分析上机第四次作业
实验项目：
共轭梯度法求解对称正定的线性方程组
实验内容：
用共轭梯度法求解下面方程组
（1）
迭代20次或满足
时停止计算。
（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次或满足
时停止计算。
（3）*考虑模型问题，方程为
用正方形网格离散化，若取
，得到
的线性方程组，并用共轭梯度法（CG法）求解，并对解作图。
实验要求：
迭代初值可以取
，计算到
停止．本题有精确解
，这里
表示以
为分量的向量，
表示在相应点
上取值作为分量的向量．
实验一：
（1）
编制函数子程序CGmethod。
function[x,k]=CGmethod（A,b）
n=length（A）;x=zeros（n,1）;r=b-A*x;rho=r'*r;
k=0;
whilerho>10^（-12）&k<20
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
编制主程序shiyan1_1：
clear,clc
A=[2,-1,0,0;-1,3,-1,0;0,-1,4,1;0,0,-1,5];
b=[3,-2,1,5]';
[x,k]=CGmethod（A,b）
运行结果为：
x=
1.3882
-0.2855
-0.0222
0.9367
k=
20
（2）
编制函数子程序CGmethod_1
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
编制主程序shiyan1_2:
clear,clc
n=1000;
A=hilb（n）;
b=sum（A'）';
[x,k]=CGmethod_1（A,b）
运行结果为：
x的值,均接近1，迭代次数k=32
实验二
实验目的：
用复化Simpson方法、自适应复化梯形方法和Romberg方法求数值积分。
实验内容：
计算下列定积分
（1）
（2）
（3）
实验要求：
（1）分别用复化Simpson公式、自适应复化梯形公式计算要求绝对误差限为
，输出每种方法所需的节点数和积分近似值，对于自适应方法，显示实际计算节点上离散函数值的分布图；
（2）分析比较计算结果。
2、实验目的：
高斯数值积分方法用于积分方程求解。
实验内容：
线性的积分方程的数值求解，可以被转化为线性代数方程组的求解问题。
而线性代数方程组所含未知数的个数，与用来离散积分的数值方法的节点个数相同。
在节点数相同的前提下，高斯数值积分方法有较高的代数精度，用它通常会得到较好的结果。
对第二类Fredholm积分方程
首先将积分区间[a,b]等分成n份，在每个子区间上离散方程中的积分就得到线性代数方程组。
实验要求：
分别使用如下方法，离散积分方程中的积分
1.复化梯形方法；2.复化辛甫森方法；3.复化高斯方法。
求解如下的积分方程
，方程的准确解为
，
并比较各算法的优劣。
实验二
1、复化Simpson方法）
输入积分区间下限0
输入积分区间上限2
输入等分份数20
输入被积函数（以x为自变量）x^6/10-x^2+x
S=
1.1619
输入积分区间下限0
输入积分区间上限1
输入等分份数20
输入被积函数（以x为自变量）x*sqrt（x）
S=
0.4000
输入积分区间下限5
输入积分区间上限200
输入等分份数20
输入被积函数（以x为自变量）1/sqrt（x）
S=
23.8218
2、自动变步长Simpson方法
函数1：
输入积分区间下限0
输入积分区间上限2
输入为课本的第几个函数（第一个这输入1）：
1
S=1.619（过程省略）
i=
19
函数2：
输入积分区间下限0
输入积分区间上限1
输入为课本的第几个函数（第一个这输入1）：
2
S=0.4（过程省略）
i=
17
函数3：
输入积分区间下限5
输入积分区间上限200
输入为课本的第几个函数（第一个这输入1）：
3
S=23.8121（过程省略）
i=
111
编制程序如下：
Clear，clc
symsx
a=input（'输入积分区间下限'）;
b=input（'输入积分区间上限'）;
n=input（'输入等分份数'）;
ff=input（'输入被积函数（以x为自变量）'）;
h=（b-a）/n;
f=inline（ff,'x'）;
sum1=0;
sum2=0;
fori=0:
n-1
sum1=sum1+f（a+i*h+0.5*h）;
end
fori=1:
n-1
sum2=sum2+f（a+i*h）;
end
fori=0:
2*n
X（i+1,1）=f（（b-a）*i/（n*2）+a）;
end
S=h/6*（f（a）+4*sum1+2*sum2+f（b））
functionS=zdsps（n）
a=0;
b=1;
h=（b-a）/4;
f=inline（'x^（3/2）','x'）;
sum1=0;
sum2=0;
fori=0:
n-1
sum1=sum1+f（a+i*h+0.5*h）;
end
fori=1:
n-1
sum2=sum2+f（a+i*h）;
end
fori=0:
2*n
x（i+1,1）=f（（b-a）*i/（n*2）+a）;
end
S=h/6*（f（a）+4*sum1+2*sum2+f（b））;
end
functionS=zpsgs（a,b,n,ff）
h=（b-a）/n;
sum1=0;
sum2=0;
sum3=0;
sum4=0;
ifff==1
f=inline（'x^6/10-x^2+x','x'）;
end
ifff==2
f=inline（'x^（3/2）','x'）;
end
ifff==3
f=inline（'1/sqrt（x）','x'）;
end
fori=0:
n-1
sum1=sum1+f（a+i*h+0.25*h）;
sum2=sum2+f（a+i*h+0.75*h）;
sum4=sum4+f（a+i*h+0.5*h）;
end
fori=1:
n-1
sum3=sum3+f（a+i*h）;
end
fori=0:
4*n
x（i+1,1）=f（（b-a）*i/（n*4）+a）;
end
S=h/（6*2）*（f（a）+4*sum1+4*sum2+2*（sum3+sum4）+f（b））;
end
clear，clc
a=input（'输入积分区间下限'）;
b=input（'输入积分区间上限'）;
ff=input（'输入为课本的第几个函数（第一个这输入1）：
'）;
fori=2:
300
S（i）=zpsgs（a,b,（i）,ff）;
S（i+1）=zpsgs（a,b,（i+1）,ff）;
ifabs（S（i+1）-S（i））<0.5*10^（-7）
break
end
end
S%所求积分值
i%所分份数
实验三
1、对常微分方程初值问题
分别使用Euler显示方法、改进的Euler方法和经典RK法和四阶Adams预测-校正算法，求解常微分方程数值解，并与其精确解
进行作图比较。
其中多步法需要的初值由经典RK法提供。
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进行试算，观察并记录计算的结果有什么特点是否发现什么不同的现象
3、
定义函数子程序为：
functionz=f（x,y）
z=-y+2*cos（x）;
return
主程序为：
clear,clc
b=pi;
a=0;
n=100;
y
（1）=1;
h=（b-a）/n;
x=a:
h:
b;
fori=1:
n
y（i+1）=y（i）+h*f（x（i）,y（i））;
end
t1=plot（x,y,'r-'）
holdon
fori=1:
n
K1=f（x（i）,y（i））;
K2=f（x（i+1）,y（i）+h*K1）;
y（i+1）=y（i）+h*（K1+K2）/2;
end
t2=plot（x,y,'b+'）
fori=1:
n
K1=f（x（i）,y（i））;
K2=f（x（i）+0.5*h,y（i）+0.5*h*K1）;
K3=f（x（i）+0.5*h,y（i）+0.5*h*K2）;
K4=f（x（i）,y（i）+h*K3）;
y（i+1）=y（i）+h*（K1+2*K2+2*K3+K4）/6;
end
t3=plot（x,y,'ko'）
fori=1:
3
K1=f（x（i）,y（i））;
K2=f（x（i）+0.5*h,y（i）+0.5*h*K1）;
K3=f（x（i）+0.5*h,y（i）+0.5*h*K2）;
K4=f（x（i）,y（i）+h*K3）;
y（i+1）=y（i）+h*（K1+2*K2+2*K3+K4）/6;
end
fori=4:
n
z=y（i）+h/24*（55*f（x（i）,y（i））-59*f（x（i-1）,y（i-1））...
+37*f（x（i-2）,y（i-2））-9*f（x（i-3）,y（i-3）））;
y（i+1）=y（i）+h/24*（9*f（x（i+1）,z）+19*f（x（i）,y（i））...
-5*f（x（i-1）,y（i-1））+f（x（i-2）,y（i-2）））;
end
t4=plot（x,y,'g*'）
t5=ezplot（'sin（x）+cos（x）',[0,pi]）
xlabel（'x轴','FontWeight','bold'）
ylabel（'y轴','FontWeight','bold'）
legend（[t1,t2,t3,t4,t5],'向前Euler法','¸改进的Euler法','经典四阶Runge-Kutta法','四阶Adams公式','精确解'）
原图为：
局部放大图为：
由图可得：
四阶Adams公式及改进的欧拉法将为精确。
升级会员 