实验一
实验项目：共轭梯度法求解对称正定的线性方程组 实验容：用共轭梯度法求解下面方程组
(1) 123421003131020141100155x x x x -?????? ? ? ?--- ? ? ?=
? ? ?-- ? ? ?-????
?? 迭代20次或满足()
(1)
1110k k x
x --∞
-<时停止计算。
编制程序：储存m 文件
function [x,k]=CGmethod(A,b)
n=length(A);x=2*ones(n,1);r=b-A*x;rho=r'*r; k=0;
while rho>10^(-11) & k<1000 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
运行程序：
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.1176 0.2353 0.5882 1.1177
(2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵，
b[1]=3, b[n]=3, b[i]=2，i=2,3,…,n-1
迭代10000次或满足()()710k k r b Ax -=-≤时停止计算。
编制程序：储存m 文件
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
运行程序：
clear,clc n=1000; A=hilb(n); b=sum(A')'; [x,k]=CGmethod(A,b)
实验二
1、 实验目的：用复化Simpson 方法、自适应复化梯形方法和Romberg 方法求数值积分。
实验容：计算下列定积分
(1) dx x x x ????
? ??+-2
02
610 (2)?10dx x x (3) ?20051dx x 实验要求：
（1）分别用复化Simpson 公式、自适应复化梯形公式计算要求绝对误差限为7102
1
-?=
ε，输出每种方法所需的节点数和积分近似值，对于自适应方法，显示实际计算节点上离散函数值的分布图；
（2）分析比较计算结果。
程序：
syms x
f=x^6/10-x^2+x %定义函数f （x ） n=input('输入所求导数阶数:')
f2=diff(f,x,n) %求f(x)的n 阶导数
（1）复化梯形 clc
syms x%定义自变量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=0.5*10^(-7) %精度要求值
a=0 %积分下限
b=2 %积分上限
x1=fminbnd(f3,1,2) %求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值
for n=2:1000000 %求等分数n
Rn=-(b-a)/12*((b-a)/n)^2*f2(x1) %计算余项
if abs(Rn) break% 符合要求时结束 end end h=(b-a)/n %求h Tn1=0 for k=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 syms x%定义自变量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值 for n=2:1000000 %求等分数n Rn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1) %计算余项 if abs(Rn) break% 符合要求时结束 end h=(b-a)/n %求h Sn1=0 Sn2=0 for k=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) 结果： dx x x x ???? ? ??+-2 02610 用复化梯形算法计算的结果 Tn= 1.6674 等分数 n= 24764 已知值与计算值的误差 R= -6.3976 用Simpson 公式计算的结果 Sn= 1.0434 等分数 n= 76 已知值与计算值的误差 R= -6.0216 ?1 dx x x 用复化梯形算法计算的结果 Tn= 0.8985 等分数 n= 1119 已知值与计算值的误差 R= -6.1665 用Simpson 公式计算的结果 Sn= 0.9406 已知值与计算值的误差 R= -6.1245 ? 200 5 1dx x 用复化梯形算法计算的结果 Tn= 23.2353 等分数 n= 1000000 已知值与计算值的误差 R= 16.1704 用Simpson 公式计算的结果 Sn= 23.2617 等分数 n= 10647 已知值与计算值的误差 R= 16.1969 分析：在处理问题时，复化Simpson 要比复化梯度计算速度要快很多。 2、实验目的：高斯数值积分方法用于积分方程求解。 实验容：线性的积分方程的数值求解，可以被转化为线性代数方程组的求解问题。而线性代数方程组所含未知数的个数，与用来离散积分的数值方法的节点个数相同。在节点数相同的前提下，高斯数值积分方法有较高的代数精度，用它通常会得到较好的结果。对第二类Fredholm 积分方程 b t a t f ds s y s t k t y b a ≤≤+=?),()(),()( 首先将积分区间[a,b]等分成n 份，在每个子区间上离散方程中的积分就得到线性代数方程组。 实验要求：分别使用如下方法，离散积分方程中的积分 1.复化梯形方法； 2.复化辛甫森方法； 3.复化高斯方法。求解如下的积分方程 t t e ds s y e e t y --=?10 )(12)(，方程的准确解为t e ， 并比较各算法的优劣。 程序结果：当迭代次数 复化梯形方法的平均误差err=0.247 复化辛甫森方法的平均误差err=0.00116 复化高斯方法的平均误差err=0.0008 复化梯形方法的平均误差err=0.00116 复化辛甫森方法和复化高斯方法的平均误差err=0 可以看出，复化高斯方法得到的结果精度最高，复化辛普森方法比复化梯形方法的精度要高，程序： clear,clc a=0;b=1; n=5; figure fun1(a,b,n); hold on a=0;b=1; n=5; figure fun2(a,b,n); hold on figure fun3(a,b,n); 编制m函数： function y=Kfun(t,s) y=2/(exp(1)-1)*exp(t); function y=ffun(t) y=-exp(t); function y=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; for i=1:n [x(2*i-1:2*i),w(2*i-1:2*i)]=GaussLegendre(x1(i),x1(i+1),2); end x=a:h/2:b; w=x; w(1)=h/6; w(2*n+1)=h/6; for i=1:n w(2*i)=2/3*h; if i w(2*i+1)=h/3; end end function [x,w]=fhtrapz(a,b,n) h=(b-a)/n; x=a:h:b; for i=1:n+1 if i==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); [V L]=eig(CM); [x ind]=sort(diag(L)); V=V(:,ind)'; w=2*V(:,1).^2; x=x*(b-a)/2+(b+a)/2; w=w*(b-a)/2; function fun1(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); for i=1:n B(i)=ffun(x1(i)); for j=1:n A(i,j)=w1(j)*Kfun(x1(i),x1(j)); end end A=eye(n)-A; y1=(A\B)'; yN=x; for i=1:n1+1 yN(i)=ffun(x(i)); for j=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('复化值','真实值'); function fun2(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); for i=1:n B(i)=ffun(x1(i)); end for i=1:n for j=1:n A(i,j)=w1(j)*Kfun(x1(i),x1(j)); end end A=eye(n)-A; y1=(A\B)'; yN=x; for i=1:n1+1 yN(i)=ffun(x(i)); for j=1:n yN(i)=yN(i)+w1(j)*Kfun(x(i),x1(j))*y1(j); end fprintf('数值解:\n') yN fprintf('精确解:\n') y0 plot(x,yN,'x',x,y0,'.'); h=legend('复化值','真实值'); function fun3(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); for i=1:n B(i)=ffun(x1(i)); end for i=1:n for j=1:n A(i,j)=w1(j)*Kfun(x1(i),x1(j)); end end A=eye(n)-A; y1=(A\B)'; yN=x; for i=1:n1+1 yN(i)=ffun(x(i)); for j=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、对常微分方程初值问题 2cos (0)1 y y x y '=-+?? =? (0)x π<< 分别使用Euler 显示方法、改进的Euler 方法和经典RK 法和四阶Adams 预测-校正算法，求解常微分方程数值解，并与其精确解cos sin y x x =+进行作图比较。其中多步法需要的初值由经典RK 法提供。 程序：子程序 function E=Euler(fun,x0,y0,h,N) x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=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') hold on function E=Euler_modify(fun,x0,y0,h,N) x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0;y(1)=y0; for n=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; for n=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); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4); end 运行以下程序 clc clear Euler('fun',0,1,0.1,100) Euler_modify('fun',0,1,0.1,100) [x,y]=Rk_N4('fun',0,1,0.1,100) plot(x,y,'-*r') hold on x1=0:0.1:pi; y1=cos(x1)+sin(x1) plot(x1,y1,'-.b') title('误差分析'); xlabel('x 轴'); ylabel('y 轴'); legend('Euler','Euler 改进法','R_K 法','精确'); axis([0 pi -1.5 1.5]); grid on 画出图形进行比较： 2、实验目的：Lorenz 问题与混沌 实验容：考虑著名的Lorenz 方程 ?????????-=--=-=)1.8()(bz xy dt dz xz y rx dt dy x y s dt dx 其中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函数程序 function ydot=blzeq(t,y) global SIGMA RHO BETA A=[-BETA 0 y(2);0 -SIGMA SIGMA;-y(2) RHO -1]; ydot=A*y; 主程序： clf clc global SIGMA RHO BETA SIGMA=10.; RHO=28.; BETA=8/3; axis([0 50 -30 30 -30 30]) view(3) hold on title('Lorenz Attractor') y0=[50,50,50]; tfinal=100; [t,y]=ode23('blzeq',[0 tfinal],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也会有发生类似的情况。
相关主题
文本预览