数值分析上机第四次作业.doc

数值分析上机第四次作业
实验项目:共轭梯度法求解对称正定得线性方程组
实验内容:用共轭梯度法求解下面方程组
(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;
while rho>10^(-12) & k<20
k=k+1;
if k==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;
while norm(r1,1)>=10^(-7)&k<10^4
k=k+1;
if k==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
syms x
a=input('输入积分区间下限');
b=input('输入积分区间上限');
n=input('输入等分份数');
ff=input('输入被积函数(以x为自变量)');
h=(b-a)/n;
f=inline(ff,'x');
sum1=0;
sum2=0;
for i=0:n-1
sum1=sum1+f(a+i*h+0、5*h);
end
for i=1:n-1
sum2=sum2+f(a+i*h);
end
for i=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))
function S = zdsps( n )
a=0;
b=1;
h=(b-a)/4;
f=inline('x^(3/2)','x');
sum1=0;
sum2=0;
for i=0:n-1
sum1=sum1+f(a+i*h+0、5*h);
end
for i=1:n-1
sum2=sum2+f(a+i*h);
end
for i=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
function S = zpsgs(a,b,n,ff )
h=(b-a)/n;
sum1=0;
sum2=0;
sum3=0;
sum4=0;
if ff==1
f=inline('x^6/10-x^2+x','x');
end
if ff==2
f=inline('x^(3/2)','x');
end
if ff==3
f=inline('1/sqrt(x)','x');
end
for i=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
for i=1:n-1
sum3=sum3+f(a+i*h);
end
for i=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):');
for i=2:300
S(i)=zpsgs(a,b,(i),ff);
S(i+1)=zpsgs(a,b,(i+1),ff);
if abs(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、
定义函数子程序为:
function z=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;
for i=1:n
y(i+1)=y(i)+h*f(x(i),y(i));
end
t1=plot(x,y,'r-')
hold on
for i=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+')
for i=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')
for i=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
for i=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公式及改进得欧拉法将为精确。