上海电力学院
数值计算方法上机实习
报告
院 系：电气工程学院 专业年级：电力系统及其自动化 学生姓名： 杨 建 学号：ys1310301074
指导教师： 黄建雄
2013年12月23日
数值计算方法上机实习题
x n
dx ， 1． 设I n =⎰05+x
1
（1） 由递推公式I n =-5I n -1+
1
，从I 0的几个近似值出发，计算I 20； n
11
I n +，计算I 0； 55n
（2） 粗糙估计I 20，用I n -1=-
（3） 分析结果的可靠性及产生此现象的原因(重点分析原因) 。
%上机习题1 %（1） I=0.182; for n=1:20
I=(-5)*I+1/n; end
fprintf('I20 的值是 %e\n',I); %(2)
I=0.0087; for n=1:20
I=(-1/5)*I+1/(5*n); end
fprintf('I0 的值是 %f\n',I);
3） 现象产生的原因：
∗∗∗∗
假设Sn的真值为Sn，误差为εn，即εn=Sn−Sn 对于真值，我们也有关系式Sn+5Sn−1=n综
1
合两个递推公式有εn=−5εn−1 这就意味着哪怕开始只有一点点误差，只要n 足够大，按照这种每计算一步误差增长5倍的方式，所得结果总是不可信的。因此整个算法是不稳定的。 对于第二种方法 εn=−εn−1 误差会以每计算一步缩小到1/5的方式进行，所以以这种方
51
法计算的结果是可靠的，整个算法是稳定的。
x
2. 求方程e +10x -2=0的近似根，要求x k +1-x k
-4
（1） 在[0，1]上用二分法； （2） 取初值x 0=0，并用迭代x k +1（3） 加速迭代的结果；
（4） 取初值x 0=0，并用牛顿迭代法； （5） 分析绝对误差。
2-e x =；
10
%（1）在[0，1]上用二分法； a=0;b=1.0;ci=0;
while abs(b-a)>5*1e-4 c=(b+a)/2; ci=ci+1;
ifexp(c)+10*c-2>0 b=c; else a=c; end end
fprintf('二分法求得 c 的值是 %f\t',c); fprintf('迭代次数是 %d\n',ci); %（2）不动点迭代法， x=0; a=1; ci=0;
while abs(x-a)>5*1e-4 a=x;
x=(2-exp(x))/10; ci=ci+1; end
fprintf('不动点迭代求得 x 的值是 %f\t',x); fprintf('迭代次数是 %d\n',ci); %（3）艾特肯加速迭代 x=0;
a=0;b=1;ci=0;
while abs(b-a)>5*1e-4 a=x;
y=exp(x)+10*x-2; z=exp(y)+10*y-2;
x=x-(y-x)^2/(z-2*y+x); b=x; ci=ci+1; end
fprintf('艾特肯加速迭代求得 x 的值是 %f\t',x); fprintf('迭代次数是 %d\n',ci); %（4）用牛顿迭代法 x=0;
a=0;b=1;ci=0;
while abs(b-a)>5*1e-4 a=x;
x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x; ci=ci+1;
end
fprintf('牛顿迭代法求得 x 的值是 %f\t',x); fprintf('迭代次数是 %d\n',ci); %（5）分析绝对误差
solve('exp(x)+10*x-2=0'); fprintf('x的真值是%.8f',x);
运行结果：
二分法求得 c 的值是 0.090332 迭代次数是 11 绝对为误差为0.000193 不动点迭代求得 x 的值是 0.090513 迭代次数是 4 绝对为误差为0.000012 艾特肯加速迭代求得 x 的值是 0.099488 迭代次数是 3 绝对为误差为0.008963 牛顿迭代法求得 x 的值是 0.090525 迭代次数是 2 绝对为误差为0.000000 x 的真值是0.090525>>
分析可知加速迭代绝对误差最大，二分法次之，牛顿法和不动点法迭代效果最好。
3．钢水包使用次数多以后，钢包的容积增大，数据如下：
试从中找出使用次数和容积之间的关系，计算均方差。（注：增速减少，用何种模型） 解：设y=f(x)具有指数形式
b
x
1
y =ae （a>0,b
x
记A=lna，B=b，并引入新变量z=lny，t=1/x。引入新变量后的数据表如下
t=[0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625];
z=[1.8594 2.1041 2.2597 2.2513 2.2721 2.3026 2.2956 2.3016 2.3504 2.3599 2.3609 2.3795 2.3609 2.3888 2.3758]; polyfit(t,z,1)
结果：
ans =-1.1107 2.4578
由此可得 A=2.4578，B=-1.1107，a =e A =11. 6791，b=B=-1.1107 方程即为y =11. 6791e
-1. 1107
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(x)=11.6791*exp( -1.1107./x); c=0;
fori=1:15 a=y(i); b=x(i);
c=c+(a-f(b))^2; end
averge=c/15
结果：averge =0.0594
0⎫⎛4-10-10⎛0⎫
⎪ ⎪ -14-10-10⎪ 5⎪ 0-14-10-1⎪ -2⎪
⎪，b = ⎪，A x =b 4．设A =
-10-14-10⎪ 5⎪ 0-10-14-1⎪ -2⎪ ⎪ ⎪ 0⎪ 6⎪0-10-14⎭⎝⎝⎭
分析下列迭代法的收敛性，并求x k +1-x k
a) JACOBI 迭代；
b) GAUSS-SEIDEL 迭代； c) SOR 迭代（ω=1. 334, 1. 95,
a) %JACOBI迭代
function y=jacobi(A,b,x0) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); B=D\(L+U); f=D\b;
y=B*x0+f;n=1;
while norm (y-x0)>1e-4
。 0. 95）
2
≤10-4的近似解及相应的迭代次数。
x0=y;
y=B*x0+f;n=n+1; end y n
b) %高斯-赛德尔迭代 function y=gs(A,b,x0) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); G=(D-L)\U; f=(D-L)\b; y=G*x0+f;n=1;
while norm(y-x0)>10^(-4) x0=y;
y=G*x0+f;n=n+1; end y n
c)%SOR迭代(omega=1.334,1.95,0.95) function y=sor(A,b,w,x0) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1);
lw=(D-w*L)\((1-w)*D+w*U); f=(D-w*L)\b*w; y=lw*x0+f;n=1;
while norm(y-x0)>10^(-4) x0=y;
y=lw*x0+f;n=n+1; end y n
A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6]'; x0=[0 0 0 0 0 0]';
fprintf('雅可比迭代的结果及迭代次数是：\n'); jacobi(A,b,x0);
fprintf('高斯赛德尔迭代的结果及迭代次数是：\n'); gs(A,b,x0);
c=[1.334 1.95 0.95]; fori=1:3
w=c(i);
fprintf('omega的值为:%.3f 时，SOR 迭代的结果及迭代次数是：\n',w) sor(A,b,w,x0); end
运行结果：
雅可比迭代的结果及迭代次数是： y = 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 n =28
高斯赛德尔迭代的结果及迭代次数是： y =1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 n =15
omega 的值为:1.334 时，SOR 迭代的结果及迭代次数是： y =1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 n = 13
omega 的值为:1.950 时，SOR 迭代的结果及迭代次数是： y =1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 n =241
omega 的值为:0.950 时，SOR 迭代的结果及迭代次数是： y =1.0000 2.0000 1.0000 2.0000 1.0000
2.0000 n = 17
⎛631⎫ ⎪
5．用逆幂迭代法求A = 321⎪最接近于11的特征值和特征向量，准确到10-3。
111⎪⎝⎭
程序：
function [mt,my]=maxtr(A,p,ep) n=length(A); B=A-p*eye(n); v0=ones(n,1); k=1; v=B*v0;
while abs(norm(v,inf)-norm(v0,inf))>ep %norm(v-v0)>ep k=k+1; q=v;
u=v/norm(v,inf) v=B*u; v0=q; end
mt=1/norm(v,inf)+p my=u
A=[6 3 1;3 2 1;1 1 1;]; maxtr(A,11,0.001); 运行结果： mt =11.0919 my =0.3845 -1.0000 0.7306
6．用经典R-K 方法求解初值问题
'=-2y 1+y 2+2sin x ⎧y 1⎧y 1(0) =2
（1）⎨，x ∈[0, 10]， ⎨；
'y =y -2y +2cos x -2sin x y (0) =312⎩2⎩2
（2）⎨
'=-2y 1+y 2+2sin x ⎧y 1⎧y 1(0) =2
，x ∈[0, 10]， ⎨。
'=998y 1-999y 2+999cos x -999sin x ⎩y 2⎩y 2(0) =3
⎧y 1(x ) =2e -x +sin x
和精确解⎨比较，分析结论。 -x
⎩y 2(x ) =2e +cos x
（1）程序：
functionydot=lorenzeq(x,y)
ydot=[-2*y(1)+y(2)+2*sin(x);y(1)-2*y(2)+2*cos(x)-2*sin(x)]
以文件名lorenzeq.m 保存。
主窗口输入：[x,y]=ode45('lorenzeq',[0:10],[2;3]) 运行结果为： x = 0 1 2 3 4 5 6 7 8 9 10
y =2.0000 3.0000 1.5775 1.2758 1.1802 -0.1457 0.2406 -0.8903 -0.7202 -0.6170 -0.9454 0.2971 -0.2745 0.9652 0.6589 0.7557 0.9901 -0.1449 0.4124 -0.9109 -0.5440 -0.8389 （2）：程序：
functionydot=lorenzeq1(x,y)
ydot=[-2*y(1)+y(2)+2*sin(x);998*y(1)-999*y(2)+999*cos(x)-999*sin(x)]; 以文件名lorenzeq1.m 保存。 程序：x=0:10;
y1=2*exp(-x)+sin(x); y2=2*exp(-x)+cos(x);
[x,y]=ode45('lorenzeq1',[0:10],[2;3]);
fprintf(' x y(1) y1 y(2) y2\n') for j=1:length(y)
fprintf('%4d %.4f %.4f %.4f %.4f\n',x(j),y(j,1),y1(j),y(j,2),y2(j)) end
运行结果为:
x y(1) y1 y(2) y2 0 2.0000 2.0000 3.0000 3.0000 1 1.5772 1.5772 1.2759 1.2761 2 1.1800 1.1800 -0.1455 -0.1455 3 0.2407 0.2407 -0.8904 -0.8904 4 -0.7202 -0.7202 -0.6169 -0.6170
5 -0.9454 -0.9454 0.2972 0.2971 6 -0.2745 -0.2745 0.9648 0.9651 7 0.6588 0.6588 0.7554 0.7557 8 0.9900 0.9900 -0.1448 -0.1448 9 0.4124 0.4124 -0.9106 -0.9109 10 -0.5439 -0.5439 -0.8389 -0.8390 结论：R-K 方法求解的结果精度较高。
7．用有限差分法求解边值问题（h=0.1）：
⎧y ''-(1+x 2) y =0
. ⎨
⎩y (-1) =y (1) =1
程序为： h=0.1;
n=(1-(-1))/h+1; x(1)=-1;x(n-1)=1; y(1)=1;y(n-1)=1; fori=1:n-1
x(i)=x(1)+(i-1)*h; q(i)=(1+x(i)^2); B(i)=2/(h^2)+q(i); end
fori=1:n-2 C(i)=-1/(h^2); end
H=diag(B)+diag(C,1)+diag(C,-1); g(1)=0+1/(h^2); g(n-1)=0+1/(h^2); fori=2:n-2 g(i)=0; end y=H\g'
运行结果为： y = 0.9027 0.8235 0.7592 0.7074 0.6661 0.6338 0.6095 0.5922 0.5814 0.5767 0.5778
0.5846
0.5974
0.6163
0.6420
0.6752
0.7167
0.7680
0.8308
0.9072
8．P100，例4.4.6用函数y=asinbx拟合数据
解：使用非线性拟合函数nlinfit ，编写MATLAB 程序如下：
function [err,a,b]=nlfitb(x,y)
if nargin
x=[1:8]'/10;
y=[0.6 1.1 1.6 1.8 2.0 1.9 1.7 1.3]';
end
beta0=[1 1]';
beta=nlinfit(x,y,@mymodel,beta0);
fprintf('The nonlinear least square fitting y=a*sin(b*x) for data\n\n'); fprintf('%6.1f',x);
fprintf('\n');
fprintf('%6.1f',y);
fprintf('\n\n is\n\t y=%7.4f*sin(%7.4f*x)\n\n',beta);
z=linspace(x(1),x(end),100);
plot(x,y,'ro' ,z,beta(1)*sin(beta(2)*z),'b-.' )
grid on ;
function yb=mymodel(beta,xb)
yb=beta(1)*sin(beta(2)*xb);
拟合结果为：
The nonlinear least square fitting y=a*sin(b*x) for data
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
0.6 1.1 1.6 1.8 2.0 1.9 1.7 1.3
is
y= 1.9750*sin( 3.0249*x)
拟合曲线为：
9．拟合形如f （x ）≈（a+bx）/(1+cx)的函数的一种快速方法是将最小二乘法用于下列问题：f （x ）(1+cx)≈（a+bx），试用这一方法拟合表4-4给出的中国人口数据。
函数是1，x ，-xf(x)而数据为（xi,f(xi)）的最小二乘拟合问题。编程如下： x=[1953 1964 1982 1900 2000]';
y=[5.82 6.95 10.08 11.34 12.66]';
A=[ones(5,1) x -x.*y];
Z=A\y;
a=Z(1)
b=Z(2)
c=Z(3)
v=linspace(1953,2000,100);
plot(x,y,'b-+',v,(a+b*v)./(1+c*v),'k-');
求得的解为a =11.5250b =-0.0059c =-5.0979e-04
10. 已知20世纪美国人口的统计数字如表（单位：百万）：
试拟合美国人口20世纪的增长率。
t=1900:10:1990;
x=[76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4];
% 2-point
p=diff(x)./diff(t);
v=p./x(1:end-1);
fprintf('%4d %6.3f %6.3f\n',[t(1:end-1);p;v]);
fprintf('\n\n\n');
% 3-point
p=(-3*x(1:end-2)+4*x(2:end-1)-x(3:end))/20;
v=p./x(1:end-2);
p2=(x(end-3:end-2)-4*x(end-2:end-1)+3*x(end-1:end))/20;
v2=p2./x(end-1:end);
p=[p p2];
v=[v v2];
fprintf('%4d %6.3f %6.3f\n',[t;p;v]);
二点公式结果：
年份 1900 1910 1920 1930 1940 1950 1960 1970 1980
年增长 1.60 1.45 1.67 0.85 1.90 2.86 2.47 2.25 2.49
年增长率 0.021 0.016 0.016 0.007 0.014 0.019 0.014 0.011 0.011 三点公式
年份 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990
年增长 1.675 1.34 2.08 0.325 1.42 3.055 2.58 2.13 2.14 2.61
年增长率 0.022 0.015 0.020 0.003 0.011 0.020 0.014 0.010 0.009 0.010
上海电力学院
数值计算方法上机实习
报告
院 系：电气工程学院 专业年级：电力系统及其自动化 学生姓名： 杨 建 学号：ys1310301074
指导教师： 黄建雄
2013年12月23日
数值计算方法上机实习题
x n
dx ， 1． 设I n =⎰05+x
1
（1） 由递推公式I n =-5I n -1+
1
，从I 0的几个近似值出发，计算I 20； n
11
I n +，计算I 0； 55n
（2） 粗糙估计I 20，用I n -1=-
（3） 分析结果的可靠性及产生此现象的原因(重点分析原因) 。
%上机习题1 %（1） I=0.182; for n=1:20
I=(-5)*I+1/n; end
fprintf('I20 的值是 %e\n',I); %(2)
I=0.0087; for n=1:20
I=(-1/5)*I+1/(5*n); end
fprintf('I0 的值是 %f\n',I);
3） 现象产生的原因：
∗∗∗∗
假设Sn的真值为Sn，误差为εn，即εn=Sn−Sn 对于真值，我们也有关系式Sn+5Sn−1=n综
1
合两个递推公式有εn=−5εn−1 这就意味着哪怕开始只有一点点误差，只要n 足够大，按照这种每计算一步误差增长5倍的方式，所得结果总是不可信的。因此整个算法是不稳定的。 对于第二种方法 εn=−εn−1 误差会以每计算一步缩小到1/5的方式进行，所以以这种方
51
法计算的结果是可靠的，整个算法是稳定的。
x
2. 求方程e +10x -2=0的近似根，要求x k +1-x k
-4
（1） 在[0，1]上用二分法； （2） 取初值x 0=0，并用迭代x k +1（3） 加速迭代的结果；
（4） 取初值x 0=0，并用牛顿迭代法； （5） 分析绝对误差。
2-e x =；
10
%（1）在[0，1]上用二分法； a=0;b=1.0;ci=0;
while abs(b-a)>5*1e-4 c=(b+a)/2; ci=ci+1;
ifexp(c)+10*c-2>0 b=c; else a=c; end end
fprintf('二分法求得 c 的值是 %f\t',c); fprintf('迭代次数是 %d\n',ci); %（2）不动点迭代法， x=0; a=1; ci=0;
while abs(x-a)>5*1e-4 a=x;
x=(2-exp(x))/10; ci=ci+1; end
fprintf('不动点迭代求得 x 的值是 %f\t',x); fprintf('迭代次数是 %d\n',ci); %（3）艾特肯加速迭代 x=0;
a=0;b=1;ci=0;
while abs(b-a)>5*1e-4 a=x;
y=exp(x)+10*x-2; z=exp(y)+10*y-2;
x=x-(y-x)^2/(z-2*y+x); b=x; ci=ci+1; end
fprintf('艾特肯加速迭代求得 x 的值是 %f\t',x); fprintf('迭代次数是 %d\n',ci); %（4）用牛顿迭代法 x=0;
a=0;b=1;ci=0;
while abs(b-a)>5*1e-4 a=x;
x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x; ci=ci+1;
end
fprintf('牛顿迭代法求得 x 的值是 %f\t',x); fprintf('迭代次数是 %d\n',ci); %（5）分析绝对误差
solve('exp(x)+10*x-2=0'); fprintf('x的真值是%.8f',x);
运行结果：
二分法求得 c 的值是 0.090332 迭代次数是 11 绝对为误差为0.000193 不动点迭代求得 x 的值是 0.090513 迭代次数是 4 绝对为误差为0.000012 艾特肯加速迭代求得 x 的值是 0.099488 迭代次数是 3 绝对为误差为0.008963 牛顿迭代法求得 x 的值是 0.090525 迭代次数是 2 绝对为误差为0.000000 x 的真值是0.090525>>
分析可知加速迭代绝对误差最大，二分法次之，牛顿法和不动点法迭代效果最好。
3．钢水包使用次数多以后，钢包的容积增大，数据如下：
试从中找出使用次数和容积之间的关系，计算均方差。（注：增速减少，用何种模型） 解：设y=f(x)具有指数形式
b
x
1
y =ae （a>0,b
x
记A=lna，B=b，并引入新变量z=lny，t=1/x。引入新变量后的数据表如下
t=[0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625];
z=[1.8594 2.1041 2.2597 2.2513 2.2721 2.3026 2.2956 2.3016 2.3504 2.3599 2.3609 2.3795 2.3609 2.3888 2.3758]; polyfit(t,z,1)
结果：
ans =-1.1107 2.4578
由此可得 A=2.4578，B=-1.1107，a =e A =11. 6791，b=B=-1.1107 方程即为y =11. 6791e
-1. 1107
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(x)=11.6791*exp( -1.1107./x); c=0;
fori=1:15 a=y(i); b=x(i);
c=c+(a-f(b))^2; end
averge=c/15
结果：averge =0.0594
0⎫⎛4-10-10⎛0⎫
⎪ ⎪ -14-10-10⎪ 5⎪ 0-14-10-1⎪ -2⎪
⎪，b = ⎪，A x =b 4．设A =
-10-14-10⎪ 5⎪ 0-10-14-1⎪ -2⎪ ⎪ ⎪ 0⎪ 6⎪0-10-14⎭⎝⎝⎭
分析下列迭代法的收敛性，并求x k +1-x k
a) JACOBI 迭代；
b) GAUSS-SEIDEL 迭代； c) SOR 迭代（ω=1. 334, 1. 95,
a) %JACOBI迭代
function y=jacobi(A,b,x0) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); B=D\(L+U); f=D\b;
y=B*x0+f;n=1;
while norm (y-x0)>1e-4
。 0. 95）
2
≤10-4的近似解及相应的迭代次数。
x0=y;
y=B*x0+f;n=n+1; end y n
b) %高斯-赛德尔迭代 function y=gs(A,b,x0) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); G=(D-L)\U; f=(D-L)\b; y=G*x0+f;n=1;
while norm(y-x0)>10^(-4) x0=y;
y=G*x0+f;n=n+1; end y n
c)%SOR迭代(omega=1.334,1.95,0.95) function y=sor(A,b,w,x0) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1);
lw=(D-w*L)\((1-w)*D+w*U); f=(D-w*L)\b*w; y=lw*x0+f;n=1;
while norm(y-x0)>10^(-4) x0=y;
y=lw*x0+f;n=n+1; end y n
A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6]'; x0=[0 0 0 0 0 0]';
fprintf('雅可比迭代的结果及迭代次数是：\n'); jacobi(A,b,x0);
fprintf('高斯赛德尔迭代的结果及迭代次数是：\n'); gs(A,b,x0);
c=[1.334 1.95 0.95]; fori=1:3
w=c(i);
fprintf('omega的值为:%.3f 时，SOR 迭代的结果及迭代次数是：\n',w) sor(A,b,w,x0); end
运行结果：
雅可比迭代的结果及迭代次数是： y = 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 n =28
高斯赛德尔迭代的结果及迭代次数是： y =1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 n =15
omega 的值为:1.334 时，SOR 迭代的结果及迭代次数是： y =1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 n = 13
omega 的值为:1.950 时，SOR 迭代的结果及迭代次数是： y =1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 n =241
omega 的值为:0.950 时，SOR 迭代的结果及迭代次数是： y =1.0000 2.0000 1.0000 2.0000 1.0000
2.0000 n = 17
⎛631⎫ ⎪
5．用逆幂迭代法求A = 321⎪最接近于11的特征值和特征向量，准确到10-3。
111⎪⎝⎭
程序：
function [mt,my]=maxtr(A,p,ep) n=length(A); B=A-p*eye(n); v0=ones(n,1); k=1; v=B*v0;
while abs(norm(v,inf)-norm(v0,inf))>ep %norm(v-v0)>ep k=k+1; q=v;
u=v/norm(v,inf) v=B*u; v0=q; end
mt=1/norm(v,inf)+p my=u
A=[6 3 1;3 2 1;1 1 1;]; maxtr(A,11,0.001); 运行结果： mt =11.0919 my =0.3845 -1.0000 0.7306
6．用经典R-K 方法求解初值问题
'=-2y 1+y 2+2sin x ⎧y 1⎧y 1(0) =2
（1）⎨，x ∈[0, 10]， ⎨；
'y =y -2y +2cos x -2sin x y (0) =312⎩2⎩2
（2）⎨
'=-2y 1+y 2+2sin x ⎧y 1⎧y 1(0) =2
，x ∈[0, 10]， ⎨。
'=998y 1-999y 2+999cos x -999sin x ⎩y 2⎩y 2(0) =3
⎧y 1(x ) =2e -x +sin x
和精确解⎨比较，分析结论。 -x
⎩y 2(x ) =2e +cos x
（1）程序：
functionydot=lorenzeq(x,y)
ydot=[-2*y(1)+y(2)+2*sin(x);y(1)-2*y(2)+2*cos(x)-2*sin(x)]
以文件名lorenzeq.m 保存。
主窗口输入：[x,y]=ode45('lorenzeq',[0:10],[2;3]) 运行结果为： x = 0 1 2 3 4 5 6 7 8 9 10
y =2.0000 3.0000 1.5775 1.2758 1.1802 -0.1457 0.2406 -0.8903 -0.7202 -0.6170 -0.9454 0.2971 -0.2745 0.9652 0.6589 0.7557 0.9901 -0.1449 0.4124 -0.9109 -0.5440 -0.8389 （2）：程序：
functionydot=lorenzeq1(x,y)
ydot=[-2*y(1)+y(2)+2*sin(x);998*y(1)-999*y(2)+999*cos(x)-999*sin(x)]; 以文件名lorenzeq1.m 保存。 程序：x=0:10;
y1=2*exp(-x)+sin(x); y2=2*exp(-x)+cos(x);
[x,y]=ode45('lorenzeq1',[0:10],[2;3]);
fprintf(' x y(1) y1 y(2) y2\n') for j=1:length(y)
fprintf('%4d %.4f %.4f %.4f %.4f\n',x(j),y(j,1),y1(j),y(j,2),y2(j)) end
运行结果为:
x y(1) y1 y(2) y2 0 2.0000 2.0000 3.0000 3.0000 1 1.5772 1.5772 1.2759 1.2761 2 1.1800 1.1800 -0.1455 -0.1455 3 0.2407 0.2407 -0.8904 -0.8904 4 -0.7202 -0.7202 -0.6169 -0.6170
5 -0.9454 -0.9454 0.2972 0.2971 6 -0.2745 -0.2745 0.9648 0.9651 7 0.6588 0.6588 0.7554 0.7557 8 0.9900 0.9900 -0.1448 -0.1448 9 0.4124 0.4124 -0.9106 -0.9109 10 -0.5439 -0.5439 -0.8389 -0.8390 结论：R-K 方法求解的结果精度较高。
7．用有限差分法求解边值问题（h=0.1）：
⎧y ''-(1+x 2) y =0
. ⎨
⎩y (-1) =y (1) =1
程序为： h=0.1;
n=(1-(-1))/h+1; x(1)=-1;x(n-1)=1; y(1)=1;y(n-1)=1; fori=1:n-1
x(i)=x(1)+(i-1)*h; q(i)=(1+x(i)^2); B(i)=2/(h^2)+q(i); end
fori=1:n-2 C(i)=-1/(h^2); end
H=diag(B)+diag(C,1)+diag(C,-1); g(1)=0+1/(h^2); g(n-1)=0+1/(h^2); fori=2:n-2 g(i)=0; end y=H\g'
运行结果为： y = 0.9027 0.8235 0.7592 0.7074 0.6661 0.6338 0.6095 0.5922 0.5814 0.5767 0.5778
0.5846
0.5974
0.6163
0.6420
0.6752
0.7167
0.7680
0.8308
0.9072
8．P100，例4.4.6用函数y=asinbx拟合数据
解：使用非线性拟合函数nlinfit ，编写MATLAB 程序如下：
function [err,a,b]=nlfitb(x,y)
if nargin
x=[1:8]'/10;
y=[0.6 1.1 1.6 1.8 2.0 1.9 1.7 1.3]';
end
beta0=[1 1]';
beta=nlinfit(x,y,@mymodel,beta0);
fprintf('The nonlinear least square fitting y=a*sin(b*x) for data\n\n'); fprintf('%6.1f',x);
fprintf('\n');
fprintf('%6.1f',y);
fprintf('\n\n is\n\t y=%7.4f*sin(%7.4f*x)\n\n',beta);
z=linspace(x(1),x(end),100);
plot(x,y,'ro' ,z,beta(1)*sin(beta(2)*z),'b-.' )
grid on ;
function yb=mymodel(beta,xb)
yb=beta(1)*sin(beta(2)*xb);
拟合结果为：
The nonlinear least square fitting y=a*sin(b*x) for data
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
0.6 1.1 1.6 1.8 2.0 1.9 1.7 1.3
is
y= 1.9750*sin( 3.0249*x)
拟合曲线为：
9．拟合形如f （x ）≈（a+bx）/(1+cx)的函数的一种快速方法是将最小二乘法用于下列问题：f （x ）(1+cx)≈（a+bx），试用这一方法拟合表4-4给出的中国人口数据。
函数是1，x ，-xf(x)而数据为（xi,f(xi)）的最小二乘拟合问题。编程如下： x=[1953 1964 1982 1900 2000]';
y=[5.82 6.95 10.08 11.34 12.66]';
A=[ones(5,1) x -x.*y];
Z=A\y;
a=Z(1)
b=Z(2)
c=Z(3)
v=linspace(1953,2000,100);
plot(x,y,'b-+',v,(a+b*v)./(1+c*v),'k-');
求得的解为a =11.5250b =-0.0059c =-5.0979e-04
10. 已知20世纪美国人口的统计数字如表（单位：百万）：
试拟合美国人口20世纪的增长率。
t=1900:10:1990;
x=[76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4];
% 2-point
p=diff(x)./diff(t);
v=p./x(1:end-1);
fprintf('%4d %6.3f %6.3f\n',[t(1:end-1);p;v]);
fprintf('\n\n\n');
% 3-point
p=(-3*x(1:end-2)+4*x(2:end-1)-x(3:end))/20;
v=p./x(1:end-2);
p2=(x(end-3:end-2)-4*x(end-2:end-1)+3*x(end-1:end))/20;
v2=p2./x(end-1:end);
p=[p p2];
v=[v v2];
fprintf('%4d %6.3f %6.3f\n',[t;p;v]);
二点公式结果：
年份 1900 1910 1920 1930 1940 1950 1960 1970 1980
年增长 1.60 1.45 1.67 0.85 1.90 2.86 2.47 2.25 2.49
年增长率 0.021 0.016 0.016 0.007 0.014 0.019 0.014 0.011 0.011 三点公式
年份 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990
年增长 1.675 1.34 2.08 0.325 1.42 3.055 2.58 2.13 2.14 2.61
年增长率 0.022 0.015 0.020 0.003 0.011 0.020 0.014 0.010 0.009 0.010