数值计算方法上机实习题
数值计算方法上机实习题
1 1
(2)丨20=0, 120 = 10000,用 I n 」
I n 」 一，计算丨0 ； 5 5n
答：第一个算法可得出
易知第一个算法每一步计算都把误差放大了
5倍，n 次计算后更是放大了 5n
倍，可靠性低。
第二个算法可得出
(1 )
由递推公式I n = -51心-，从
n
I o =O.1823出发，计算* ;
(3)
分析结果的可靠性及产生此现象的原因
重点分析原因。
可以看出第二个算法每一步计算就把误差缩小5倍，n次后缩小了5n倍，可靠性高。
2.求方程e^+IOx —2=0的近似根，要求 x k *—x k <^10-，并比较计算量。
(1)
在［0 , 1］上用二分法；
(1) [0, 1]上的二分法
二分法子程序：
function [root,n]=EFF3(f,x1,x2) %第二题⑴二分法
f1= subs(f,symvar(f),x1);% 函数在 x=x1 的值 f2=subs(f,symvar(f),x2);%x=x2 n=0;%步数 er=5*10A-4;% 误差 if(f1==0)
root=x1; return;
elseif(f2==0)
root=x2; return;
elseif(f1*f2>0)
disp('两端点函数值乘积大于 0！');
return;
else
while(abs(x1-x2)>er)% 循环 x3=(x1+x2)/2;
f3=subs(f,symvar(f),x3);
n=n+1; if(f3==0)
root=x3; break; elseif(f1*f3>0)
x1=x3; else
x2=x3; end end
root=(x1+x2)/2;%while 循环少一步需加上 end
、
2 _e (2) 取初值x 0 =0，并用迭代x k d :
10
(2)初值x o =O 迭代
计算根与步数程序：
迭代法子程序：
syms x;
function [root,n]=DDF(g,x0,err,max) (接下页) f=(2-exp(x))/10;(接下页)
计算根与步数程序：
fplot(@(x) exp(x)+10*x-2,[0,1]); grid on; syms x; f=exp(x)+10*x-2; [root,n]=EFF3(f,0,1);
fprintf('root=%6.8f ,n=%d \n',root,n);
计算结果显示：
root=0.09057617 ,n=11
(3)加速迭代的结果
(4)取初值X。=0，并用牛顿迭代法;
(5)分析绝对误差。
答：可以看到，在同一精度下，二分法运算了11次，题设迭代算式下运算了4次，加速迭代下运算了2次，牛顿迭代下运算了2次。因不动点迭代法和二分法都是线性收敛的，但二分法压缩系数比题设迭代方法大，收敛速度较慢。加速迭代速度是超线性收敛，牛顿法是二阶，收敛速度快。
3.钢水包使用次数多以后，钢包的容积增大，数据如下:
ax + b
试从中找出使用次数和容积之间的关系，计算均方差。(用y 拟合)
c + x
拟合曲线程序：
x=2:16;
y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59
10.60 10.8 10.6 10.9 10.76];
[f,fval]=fminsearch(@delta1,[0,0,0],optimset,x,y)
J
%fminsearch为求解多元函数的最小值函数
%f为多元函数初值x0附近的极小值点
%fval为极小值
k=2:100;
y1= (f(1).*k+f(2))./(f(3)+k); figure1=figure('Color',[1 1 1]);
gxt=plot(x,y,'xr',k,y1,'k-');
legend(原数据’,'拟合曲线','Location','northwest'); %y 为数据点连接曲线，y1为拟合曲线
title('函数y=(ax+b)/(x+c)的拟合
','FontSize',14,'FontWeight','Bold'); xlabel('次数
','FontSize',14,'FontWeight','Bold'); ylabel('容积
','FontSize',14,'FontWeight','Bold');
set(gxt,'LineWidth',1.5);
grid on;
%计算均方差
for i=1:15 构造函数子程序：
function delta=delta1(f,x,y)
a=f(1);
b=f(2);
c=f(3);
delta=0;
for k=1:15
delta=delta+((x(k)+c)*y(k)-(a*x(k)+b))A2;
end
计算结果显示：
拟合出的方程为：(x+-0.7110)y=11.2657x+-15.5024
均方差为033165089
总结：
指标选择，因题设方程I & E需［片—签普「为非线性的，要转化为线性方程故需提指标为：
y2(i)=(f(1).*x(i)+f(2))./(f(3)+x(i)); (接下页)
end j=0； for i=1:15
j=i+(y(i)-y2(i))A
2； end
jfc=sqrt(j/15);
fprintf('拟合出的方程
为:(x+%4.4f)y=%4.4fx+%4.4f \n 均方差 为:%4.8f \n',f(3),f(1),f(2),jfc);
i=i
亦2
+ l)y￡— oxf - 力]二
1=1
(cx^-i)y=ax+b f l= J-jJfcXj +1)给一叫 一 b] 其驻点
方程为：
、亠
￡ =》2[宓+Dyj-oXi- >](-%!)= 0；
i=i
计算结果显示:
12
6
10
9
8
0 10 20 30 40 50 60 70 80 90 100
次数
『4
-1 0 -1 0 0、
『0、
-1 4 -1 0 -1 0
5
0 -1 4 -1 0 -1 ，b = -2 -1 0 -1 4 -1 0
5
0 -1 0 -1 4 -1
-2 1° 0 -1 0 -1 4」
3」
(1) JACOBI 迭代;
(2) GAUSS-SEIDEL 迭代;
分析下列迭代法的收敛性，并求
X
k 1 -x
k
<10^的近似解及相应的迭代次数。
（3）SOR迭代（取=0.1:0.1:1.9，找到迭代步数最少的* ）。
逆幕迭代法子程序：
function [lamda,V ,k]=NMSF2(A,v,eps,p) [n,n]=size(A); A=A-p*eye( n); v0=v;
[tmax,tindex]=max(abs(v0)); lamd0=v0(tindex); u0=v0/lamd0; flag=0; k=0;
while(flag==0)
V=A\u0;
[tmax,tindex]=max(abs(V)); lamd1=V(tindex); u0=V/lamd1;
if (abs((lamd0)A(-1)-(lamd1)A(-1)))v=eps flag=1; end
lamd0=lamd1; k=k+1;
end (接下页) lamda=(lamd1)A(-1)+p; V=u0';
逆幕迭代计算程序：
A=[6 3 1;3 2 1;1 1 1]; v=[1 1 1]'; eps=0.001; p=11;
[lamda,U,k]=NMSF2(A,v,eps,p);
fprintf('最接近11的特征值为:%4.8f\n 特征向量： \n%4.8f\n %4.8f\n%4.8f\n',lamda,U(1),U(2),U(3));
计算结果显示：
最接近于11的特征值为：7.87313712 特征向量：
1.00000000 0.54922698 0.22545618
6.用经典R-K 方法求解初值问题
"6 3
r
5?用逆幕迭代法求 A =
3 2 1 最接近于
1 1 1
11的特征值和特征向量，准确到
10’。
'y； = _2y
”2 = yi —2y2 +2cosx—2sin x
和精确解丿y i（x） =2e」Minx比较，进行误差分析
得到结论，图形显示精确解和数值解。y2（x） = 2e? +cosx
经典四阶R-K方法
R-K程序：
f=@(x,y)
[-2*y(1)+y(2)+2*sin(x);y(1)-2*y(2)+2*cos(x)-2*sin(x)];
g=@(x,y)【-2*y(1)+y(2)+2*sin(x);998*y(1)-999*y(2)+
999*cos(x)-999*sin(x)]; (接下页)
x=0:0.1:10;
Y1= 2*exp(-x)+sin(x);
Y2=2*exp(-x)+cos(x);
x仁0:10/10000:10;
Y3=2*exp(-x1)+sin(x1);
Y4=2*exp(-x1)+cos(x1);
N=100;N 1=10000;
[a1,b1]=ode23(f,[0:0.1:10],[2 3]);
[a2,b2]=ode23(g,[0:1/1000:10],[2 3]);
[a3,b3]=ode45(f,[0:0.1:10],[2 3]);
[a4,b4]=ode45(g,[0:1/1000:10],[2 3]);
n=length(b1);
figure1=figure('Color',[1 1 1]);
dbt=plot(x,b1(:,1),'r-',x,b1(:,2),'k-',x,Y1,'g--',x,Y2,'y--'); legend('y1','y2','精确解Y1','精确解Y2');
title('四阶RK算法下的方程组的数值解&方程的精确
解','FontSize',14,'FontWeight','Bold');
ylabel( Y,'FontSize',14,'FontWeight','Bold');
xlabel('X','FontSize',14,'FontWeight','Bold');
set(dbt,'LineWidth',1.5);grid on;
Q=b1(:,1)-Y1';W=b1(:,2)-Y2';
Q1= b2(:,1)-Y3';W 仁b2(:,2)-Y4';
Q2=b3(:,1)-Y1';W2=b3(:,2)-Y2';
Q3=b4(:,1)-Y3';W3=b4(:,2)-Y4';
ER1=max(abs(Q));ER2=max(abs(W));
ER3=max(abs(Q1));ER4=max(abs(W1));
fprintf（'\n （1）中y1在ode23下求出的最大误差
为：%4.8f\n（ 1 ）中y2在ode23下求出的最大误差
为：%4.8f\n',ER1,ER2）;
fprintf（'\n （2）中y1在ode23下求出的最大误差
为：%4.8f\n（2）中y2在ode23下求出的最大误差
为：%4.8f\n',ER3,ER4）;
fprintf（'\n （1）中y1在ode45下求出的最大误差
为：%4.8f\n （ 1 ）中y2在ode45下求出的最大误差
为：%4.8f\n',ER5,ER6）;
fprintf（'\n （2）中y1在ode45下求出的最大误差
为：%4.8f\n（2）中y2在ode45下求出的最大误差
为：%4.8f\n',ER7,ER8）;
计算与显示程序：
（1 ）中y1在ode23下求出的最大误差为：
0.00169850
（1 ）中y2在ode23下求出的最大误差为：
0.00207422
（2）中y1在ode23下求出的最大误差为：
0.00000583
（2）中y2在ode23下求出的最大误差为：
0.00581329
（1 ）中y1在ode45下求出的最大误差为：
0.00052155
（1 ）中y2在ode45下求出的最大误差为：
0.00045793
（2）中y1在ode45下求出的最大误差为：
0.00000276
（2）中y2在ode45下求出的最大误差为：
0.00275331
总结：
对比2种方程在ode23,ode45命令下，与精确值的误差可以发现，在此题中，ode45的误差相对
于ode23的误差小，较适合用在此类钢性问题中。
x [0,10], 伙
0) = 2
y
2 (0) = 3
f F y i --2y1 y2 2sin x
= 998y1 -999y2 999cosx -999sinx
x [0,10],
"(0) = 2
y2(0)=3°
ER7=max(abs(Q3));ER8=max(abs(W3));
计算结果显示：
(1) Ode 23 的题1
四阶良K算法卜的方程组的数值解&方程的粘确解
四阶RK算法下的方程组的数值鞄方程的精确解
(3) Ode 45 的题1
四阶RK篦法Z的方程组的数值解&方程的秸确解
7?用有限差分法求解边值问题(h=0.1),并图形显示。
y”_(1 +X2)y =0
y(_1) = y(1)=1
有限差分法子程序：
function y=YXCFF_D(h)
n=(1-(-1))/h;
x(1)=-1;
for i=1:n-1 % n-1
x(i)=x(1)+(i-1)*h;
q(i)=(1+x(V2);
B(i)=2/(h A2)+q(i);
end
for i=1:n-2 % n-2
C(i)=-1/(hA2);
end
b(1)=0+1/(hA2);
b(n-1)=0+1/(hA2);
for i=2:n-2
b(i)=0;
end
b=b';% n-1
u(1)=b(1)/B(1);
v(1)=C(1)/B(1);
for k=2:n-2;
u(k)=(b(k)-u(k-1)*C(k-1))/(B(k)-v(k-1)*C(k-1)); v(k)=C(k)/(B(k)-v(k-
1)*C(k-1));
end
u( n-1)=(b( n-1)-u( n-2)*C( n-2))/(B( n-1)-v( n-2)*C( n-2)); V(n-1)=u(n-1); for k=n-2:-1:1
V(k)=u(k)-v(k)*V(k+1);
end
y=[1,V,1]; 计算与显示程序：
h=0.0001;
h1=0.1;
L=YXCFF_D(h);
L仁YXCFF_D(h1);
X=-1:h:1;
X仁-1:h1:1;
figure1=figure('Color',[1 1 1]);
gxt=plot(X,L,'k-',X1,L1,'r-');
title(' 有限差分法','FontSize',14,'FontWeight','Bold');
xlabel(' 输入值x','FontSize',14,'FontWeight','Bold');
ylabel('值y','FontSize',14,'FontWeight','Bold'); set(gxt,'LineWidth',1.5);
grid on;
总结：
对比两个步长下，求边值问题，步长小，越精确。
有限差分法