数值分析上机题Matlab--东南大学出版社（牛顿迭代/逐次超松弛迭代/3次样条插值/复合梯形SimpsonRomberg/四阶经典Runge-Kutta/幂法求特征向量）

第二章上机题 Newton迭代法

 

function [x,err]=Newton(f,x0,epsilon)
%用例：[x,err]=Newton('x^3/3-x',0.7,0.005)
%Input  - f  字符串公式'x^3/3-x'
%       - x0 迭代初值
%       - epsilon 是迭代精度要求
%Output – x 是最后迭代的近似结果
%       - err 是最后得到的误差
syms x
f=str2sym(f);
f(x)=f;
df(x)=diff(f(x));
phi(x)=x-f(x)/df(x);
restrain=1;
count = 0;
e = 1;
while abs(e)>epsilon
    x1=phi(x0);
    e=x1-x0;
    x0=x1;
    count = count+1;
    fprintf('已迭代%d次,', count)
    fprintf('x位：%f,', x0)
    fprintf('误差为：%f\n', e)
    if count>100
        fprintf('牛顿迭代发散\n')
        restrain=0;
        break
    end
end
if restrain==1
    fprintf('牛顿迭代收敛结束\n')
    fprintf('Newton迭代的近似解 x = %f\n',x0)
    fprintf('迭代次数count = %d\n',count)
    fprintf('误差为：%f\n', e)
end
err=vpa(e);
x=vpa(x0);
end

 

解：利用二分法在[0.7,0.8]中寻找δ，最终保留四位有效数字，得到δ=0.7625，运行结果如图1所示：

图1 二分法求解结果

②试取若干初始值，观察当x0∈-∞,-1，-1,-δ，-δ,δ，δ,1 ，1,+∞ 时Newton序列是否收敛以及收敛于哪一个根。

当x0∈-∞,-1时，Newton序列收敛，收敛于x1*=-3；

当x0∈-1,-δ时，Newton序列收敛，收敛于x1*=-3或者x3*=3；当x0取-0.81时收敛于3，而当x0取-0.8时收敛于-3；

当x0∈-δ,δ时，Newton序列收敛，收敛于x2*=0；

当x0∈δ,1时，Newton序列收敛，收敛于x1*=-3或者x3*=3；当x0取0.81时收敛于-3，而当x0取0.8时收敛于3；

当x0∈1,+∞时，Newton序列收敛，收敛于x3*=3；

 

在第二题中通过取不同的初值，即使相邻很近，但会收敛到不同的相距较远的根，我将函数曲线可视化出来，如图2所示，我的想法是，在-1的领域中，曲线斜率变化很小，会使得迭代方程的结果跳跃到别处，但由于这个函数整体是收敛的，故迭代值的跳跃并没有影响其可以得到一个结果，但这种函数我认为是不稳定的，如果在工程实践中，有这样的一个问题，输入初值是测量值，但由于各种外界因素的影响使其受到很小扰动时，得到的结果是不同的，这样是不希望发生的。所以我认为，应该尽量将方程的根求出来，当可行性较差时，也应该分析系统的灵敏性。

图2 x^3/3-x函数图像

 第三章上机题 逐次超松弛迭代法

function [x,e,count] = SOR(A,b,x,omg)
%SOR 用例
%A=[31,-13,0,0,0,-10,0,0,0;-13,35,-9,0,-11,0,0,0,0;0,-9,31,-10,0,0,0,0,0;0,0,-10,79,-30,0,0,0,-9;0,0,0,-30,57,-7,0,-5,0;0,0,0,0,-7,47,-30,0,0;0,0,0,0,0,-30,41,0,0;0,0,0,0,-5,0,0,27,-2;0,0,0,-9,0,0,0,-2,29];
% b=[-15,27,-23,0,-20,12,-7,7,10]';
% x=[0,0,0,0,0,0,0,0,0]';
% omg=0.02;
% [x,e,count] = SOR(A,b,x,omg)
%Input  - A  系数矩阵
%       - b  结果向量
%       - x  初始迭代向量
%       - omg松弛因子
%Output - x  是迭代结果
%       - e  最后最大误差
%       - count 迭代次数
 
n=length(x);
z=x;
X=[x];
e=1;
count=0;
while e>0.000005
    for i=1:n 
        summ=0;
        for j=1:n
            if i~=j
                summ=summ+A(i,j)*z(j);
            end
        end
        z(i)=(1-omg)*z(i)+omg*( b(i)-summ )/A(i,i);
    end
    X=[X,z];
    e=max(abs(z-x));
    x=z;
    count=count+1;
    if e>1000
        count=NaN;
        fprintf('超松弛因子为：%f，迭代发散\n', omg)
        break
    end
end
if flag==1
    fprintf('超松弛因子为：%f，迭代收敛，', omg)
    fprintf('迭代次数为：%d。 ', count)
end
end

 

 

运用第一问中编好的通用程序再次编程，解答程序如下：

clear,clc
Count=[];
X=[];
E=[];
A=[31,-13,0,0,0,-10,0,0,0;-13,35,-9,0,-11,0,0,0,0;0,-9,31,-10,0,0,0,0,0;0,0,-10,79,-30,0,0,0,-9;0,0,0,-30,57,-7,0,-5,0;0,0,0,0,-7,47,-30,0,0;0,0,0,0,0,-30,41,0,0;0,0,0,0,-5,0,0,27,-2;0,0,0,-9,0,0,0,-2,29];
b=[-15,27,-23,0,-20,12,-7,7,10]';
 
for i=1:99
    x=[0,0,0,0,0,0,0,0,0]';
    omg=i/50;
    [x,e,count] = SOR(A,b,x,omg);
    X=[X,x];
    E=[E,e];
    Count=[Count,count];
end
Omg=0.02:0.02:1.98;
plot(Omg,Count)
BestC=find(Count==min(Count));
fprintf('最佳松弛因子为\n')
BestOmg=Omg(BestC)

取0向量为初始迭代向量，运行结果如图3所示：

 

图3 SOR收敛过程

将结果汇总如下表：

松弛因子	迭代次数	松弛因子	迭代次数	松弛因子	迭代次数	松弛因子	迭代次数
0.02	855	0.52	42	1.02	15	1.52	22
0.04	477	0.54	40	1.04	15	1.54	24
0.06	336	0.56	38	1.06	14	1.56	25
0.08	261	0.58	37	1.08	13	1.58	26
0.10	214	0.60	35	1.10	13	1.60	28
0.12	182	0.62	34	1.12	12	1.62	31
0.14	158	0.64	32	1.14	11	1.64	34
0.16	140	0.66	31	1.16	11	1.66	36
0.18	125	0.68	30	1.18	10	1.68	39
0.20	113	0.70	29	1.20	10	1.70	44
0.22	103	0.72	27	1.22	10	1.72	50
0.24	95	0.74	26	1.24	11	1.74	55
0.26	88	0.76	25	1.26	11	1.76	63
0.28	82	0.78	24	1.28	12	1.78	76
0.30	76	0.80	24	1.30	13	1.80	89
0.32	71	0.82	23	1.32	13	1.82	110
0.34	67	0.84	22	1.34	14	1.84	144
0.36	63	0.86	21	1.36	15	1.86	207
0.38	59	0.88	20	1.38	16	1.88	361
0.40	56	0.90	19	1.40	16	1.90	1412
0.42	53	0.92	19	1.42	17	1.92	发散
0.44	51	0.94	18	1.44	17	1.94	发散
0.46	48	0.96	17	1.46	19	1.96	发散
0.48	46	0.98	17	1.48	20	1.98	发散
0.50	44	1.00	16	1.50	20

图4 松弛因子与迭代次数关系

如表格和图4所示，最佳松弛因子为1.18、1.20、1.22，迭代次数为10，解向量为：

可以看出，比ω=1时的Gauss-Seidel迭代法迭代次数16次要更优，但当ω取值较小或大时会迭代上千次，甚至发散。

第四章上机题 3次样条插值函数

function yy=cubic_spline(x,y,xx)
% 用例
% x=[0,1,2,3,4,5,6,7,8,9,10];
%y=[2.51,3.30,4.04,4.70,5.22,5.54,5.78,5.40,5.57,5.70,5.80];
% xx=[0.5:1:9.5];
% yy=cubic_spline(x,y,xx)
%Input  - x  已知数据的x坐标
%       - y  已知数据的y坐标
%       - xx 预测位置数据
%Output - yy 三次样条插值输出结果
 
%第一型，已知两端点处f(x)的1阶导数
f1(1)=0.8; f1(2)=0.2;
n=length(x); h=[];
for i=1:n-1
    h(i)=x(i+1)-x(i);
end
f=[];
for i=1:n-1
    f(i)=(y(i+1)-y(i))/h(i);
end
miu=[];
for i=1:n-2
    miu(i)=h(i)/(h(i)+h(i+1));
end
d=[];
for i=0:n-1
    if i==0
        d(i+1)=(f(i+1)-f1(1))/(h(i+1));
    elseif i==n-1
        d(i+1)=(f1(2)-f(i))/(h(i));
    else
        d(i+1)=(f(i+1)-f(i))/(h(i)+h(i+1));
    end
end
A=[];
for i=1:n
    A(i,i)=2;
    if i==1
        A(i,i+1)=1;
    elseif i==n
        A(i,i-1)=1;
    else
        A(i,i+1)=1-miu(i-1);
        A(i,i-1)=miu(i-1);
    end
end
B=6.*d';
M=A\B;  %求解M向量
 
%构造最后的插值分段函数
s2=[];
for i=1:n-1
    s2(i)=f(i)-(1/3*M(i)+1/6*M(i+1))*h(i);
end
s3=[];
for i=1:n-1
    s3(i)=1/2*M(i);
end
s4=[];
for i=1:n-1
    s4(i)=1/6/h(i)*(M(i+1)-M(i));
end
S=[y(1:n-1)',s2',s3',s4'];
a=x;
digits(5) 
for i=1:n-1
    syms x
    l1(x) = S(i,2) * (x - a(i));
    l2(x) = S(i,3) * (x - a(i)) ^ 2;
    l3(x) = S(i,4) * (x - a(i)) ^ 3;
    fprintf("-------------------------------------\n")
    Sx(x) = S(i,1) + l1 + l2  +l3;
    Sx(x) = vpa(Sx);
    fprintf('S(x)=%s \t x∈(%d,%d)\n',char(Sx),i-1,i);
end
fprintf("\n\n\n")
yy=[];
for i=1:length(xx)
    sec=xx(i)-a;
    for j=1:length(sec)
        if sec(j)<0
            flag=j-1;
            break
        end
    end
    yy(i)=S(flag,1)+S(flag,2)*(xx(i)-a(flag))+S(flag,3)*(xx(i)-a(flag))^2+S(flag,4)*(xx(i)-a(flag))^3;
end
end

 

clear,clc
x=[0,1,2,3,4,5,6,7,8,9,10];
y=[2.51,3.30,4.04,4.70,5.22,5.54,5.78,5.40,5.57,5.70,5.80];
xx=[0.5:1:9.5];
yy=cubic_spline(x,y,xx)
 
yy2 = spline(x,y,xx)
plot(xx,yy,'--',xx,yy2,'*')

 

结果展示如图5所示：

图5 3次样条插值函数

如图5所示为3次样条插值函数构造的分段函数，左图最后输出的yy为所编程序计算出的插值结果，yy2是MATLAB内嵌函数的计算结果，将两组结果绘制出来得到右图，为了方便观察，使用蓝线表示本人自己编制程序的结果，红星表示内嵌函数的结果，如果都用点，会发现两组值基本重合，说明本人所编制程序的可靠性，同时进一步研究可知，使用的内嵌函数默认是自然样条插值，属于第二型3次样条插值，将本人编制的程序进一步改进可以得到相同的结果。在编制3次样条插值函数的过程中，进一步提高了编程的逻辑构造能力。

第五章上机题

function T=trapezc(f,a,b,div)
%用例：T=trapezc@(x)sin(x)/x,0,1,32)
%Input  - f  字符串公式'sin(x)/x'
%       - a  积分左端点
%       - b  积分右端点
%       - div 确定几等分
%Output - T  输出数值积分值
 
h=(b-a)/div; %区间长度
ya=feval(f,a);
yb=feval(f,b);
%对端点值是否计算出错进行判断
if isnan(ya)==1
    ya=feval(f,a+0.001);
end
if isnan(yb)==1
    yb=feval(f,b-0.001);
end
%计算中间节点
yx=0;
for i=1:div-1
    yx=yx+feval(f,a+h*i);
end
%输出数值积分值
T=h/2*(ya+yb+2*yx);
end

function T=Simpc(f,a,b,div)
%用例：T=Simpc(@(x)sin(x)/x,0,1,32)
%Input  - f  字符串公式'sin(x)/x'
%       - a  积分左端点
%       - b  积分右端点
%       - div 确定几等分
%Output - T  输出数值积分值
 
h=(b-a)/div; %区间长度
ya=feval(f,a);
yb=feval(f,b);
%对端点值是否计算出错进行判断
if isnan(ya)==1
    ya=feval(f,a+0.001);
end
if isnan(yb)==1
    yb=feval(f,b-0.001);
end
%计算中间节点
yx1=0;
for i=1:div-1
    yx1=yx1+feval(f,a+h*i);
end
yx2=0;
for i=1:div
    yx2=yx2+feval(f,a+h*(i-0.5));
end
%输出数值积分值
T=h/6*(ya+yb+2*yx1+4*yx2);
end


function Tm=Romberg(f,a,b,k)
%用例：T=Romberg(@(x)sin(x)/x,0,1,5)
%Input  - f  字符串公式'sin(x)/x'
%       - a  积分左端点
%       - b  积分右端点
%       - k  等分2的k次方
%Output - T  输出数值积分值
 
ya=feval(f,a);
yb=feval(f,b);
%对端点值是否计算出错进行判断
if isnan(ya)==1
    ya=feval(f,a+0.001);
end
if isnan(yb)==1
    yb=feval(f,b-0.001);
end
%计算中间节点
T=[]; h=[];
T(1,1)=(b-a)/2*(ya+yb);
for i=1:k
    if i==1
        h(i)=b-a;
    else
        h(i)=h(i-1)/2;
    end
    yx=feval(f,a+0.5*h(i));
    n=2^(i-1)-1;
    for j=1:n
        yx=yx+feval(f,a+j*h(i)+0.5*h(i));
    end
    T(i+1,1)=T(i,1)/2+h(i)*yx/2;
end
n=length(T);
for m=1:n-1
    for i=m+1:n
        T(i,m+1)=((4^m)*T(i,m)-T(i-1,m))/((4^m)-1);
    end
end
if k==0
    Tm=T(1,1);
else
    Tm=T(n,m+1);
end
end

 

clear,clc
I=0.94608307036718301494;
%复合梯形公式
Tt=[1;2;4;8;16;32];
for i=1:6
    Tt(i,2)=trapezc(@(x)sin(x)/x,0,1,2^(i-1));
    Tt(i,3)=I-Tt(i,2);
end
%复合Simpson公式
Ts=[1;2;4;8;16;32];
for i=1:6
    Ts(i,2)=Simpc(@(x)sin(x)/x,0,1,2^(i-1));
    Ts(i,3)=I-Ts(i,2);
end
%Romberg积分法
Tr=[1;2;4;8;16;32];
for i=1:6
    Tr(i,2)=Romberg(@(x)sin(x)/x,0,1,i-1);
    Tr(i,3)=I-Tr(i,2);
end
fprintf('复合梯形公式\n 等分数 \t  估算值 \t   绝对误差\n');
Tr
fprintf('复合Simpson公式\n 等分数 \t  估算值 \t 绝对误差\n');
Ts
fprintf('Romberg积分法\n 等分数 \t  估算值 \t 绝对误差\n');
Tr

 

输出积分值及绝对误差如图6所示：

图6 3种数值积分方法结果与误差

第六章上机题

clear,clc
y(1)=1; x(1)=0;
h=0.1;
for i=1:31
    k1=y(i)-exp(x(i))*cos(x(i));
    k2=(y(i)+h*k1/2)-exp(x(i)+h/2)*cos(x(i)+h/2);
    k3=(y(i)+h*k2/2)-exp(x(i)+h/2)*cos(x(i)+h/2);
    x(i+1)=x(i)+h;
    k4=(y(i)+h*k3)-exp(x(i+1))*cos(x(i+1));
    y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);
end
sol=dsolve('Dy-y+exp(x)*cos(x)','y(0)=1','x');
ezplot(sol)
hold on
plot(x,y,'r')

 

用MATLAB自带的dsolve函数进行微分方程的求解，并绘制图像，与四阶经典Runge-Kutta方法求解结果进行对比，如下图所示，其中，蓝色线条为dsolve函数求解结果，红色线条为Runge-Kutta结果，结果表明龙格方法的精确性：

图7 四阶经典Runge-Kutta求解结果

 第七章上机题

clear,clc
rng(0);
A = rand(500);
v(:,1) = ones(500,1);
v(:,2)=A*v(:,1);
m(1)=max(v(:,2));
v(:,2)=v(:,2)/m(1);
n=max(abs(v(:,2)-v(:,1)));
i=3;
while n>0.00000000001
    v(:,i)=A*v(:,i-1);
    m(i-1)=max(v(:,i));
    v(:,i)=v(:,i)/m(i-1);
    n=max(v(:,i)-v(:,i-1));
    i=i+1;
end
fprintf('幂法计算的最大特征值为：%s \n',m(i-2));
fprintf('使用eig计算得到的最大特征值为：%s \n',max(eig(A)));
fprintf('最大特征值对应的特征向量为：º\n');
v(:,i-1)*m(i-2)

 

计算结果如图8所示：

图8 矩阵特征值计算结果

由结果可知，使用幂法计算得到的矩阵最大特征值与使用eig函数计算得到的最大特征值结果一致，证明了幂法程序的正确性。