华科计算方法（数值计算）上机实验题

秒速五厘米 2022-05-24 05:58 274阅读 0赞

# 《计算方法》上机实验试题 #

注：以下代码均通过matlab语言实现

1.      (25分)计算积分

![70][]

若用下列两种算法

![70 1][]

试依据积分In的性质及数值结果说明何种算法更合理。

测试代码question1.m：

format short e   
    	  
    % 算法A：由1-20开始递推  
    I_forward=zeros(21,1);  
    I_forward(1)=log(6/5); % 先计算I0，但MATLAB中下标必须从1开始  
    for n=2:21  
        I_forward(n)=-5*I_forward(n-1)+1/(n-1); % 课本上1/n要改为1/(n-1)  
    end  
    fprintf('In calculated from n=1 to n=20:\n');  
    I_forward(2:21)  
    	  
    % 算法B：由20-1开始递推  
    I_backward=zeros(21,1);  
    I_backward(21)=0.5*(1/(6*21)+1/(5*21)); % 由放缩不等式计算I20，这里下标为21 
    for n=21:-1:2  
        I_backward(n-1)=-(1/5)*I_backward(n)+1/(5*(n-1)); % 1/n同样要改为1/(n-1)  
    end  
    fprintf('In calculated from n=19 to n=0:\n');  
    I_backward(20:-1:1)

计算结果（这里每道题结果均为MATLAB截图）：

![20180522124740942][]![20180522124813813][]

左图对应算法A，从上到下为I1-I20（I0为手工计算）；右图对应算法B，从上到下为I19-I0（I20为手工计算）。

性质：In>0且单调递减。算法A从n=18开始并未出现单调递减，其每一步舍入误差都增长5倍；而算法B每一步舍入误差都减少1/5。所以算法B合理。

2.      (25分)求解方程f(x)=0有如下牛顿迭代公式

![70 2][]

(1)                 编制上述算法的通用程序，并以  ![70 3][]（ε为预定精度）作为终止迭代的准则。

(2)                 利用牛顿迭代求解方程，设预定精度ε=10-10。

![70 4][]

向量型牛顿迭代法的通用程序Newton\_iteration.m，传入的参数均可以是向量：

% newton_iteration: 向量型Newton迭代法  
    function [x1,k] = newton_iteration(f,f_de,varible,x0,eplison)  
        % dx 表示微分，f_de表示f的导数，variable表示f中的变量  
        % subs函数将f和f_de中的变量替换为x0中的值  
        dx = -subs(f_de,varible,x0)\subs(f,varible,x0);% A\b相当于inv(A)*B  
        x1=x0+dx;  
        k=1; % 迭代次数  
        while norm(x1-x0)>=eplison  
            x0=x1;  
            dx=-subs(f_de,varible,x0)\subs(f,varible,x0);  
            x1=x0+dx;  
            k=k+1;  
        end  
    end

测试代码question2.m：

format short  
    syms x y  
    f1=5*cos(x)-x*y-3;  
    f2=y^2+8*x*y-7;  
    F=[f1;f2];  
    F_de=[diff(F,'x'),diff(F,'y')]; % 这里生成4*4的雅克比矩阵  
    eplison=10e-10;  
    variable=[x;y];  
    X0=[1;1];  
    
    [x,k]=newton_iteration(F,F_de,variable,X0,eplison);  
    fprintf('after %d newton_iteration, the root is:\n', k);  
    for i=1:length(x)  
        fprintf('%1.8f\n',x(i));  
    end

计算结果：

![20180522124836936][]

由于eplison=10e-10计算时间过长，上面是将精度改成10e-4时的计算结果。

3.      (25分) 已知![70 5][]，

(1) 利用插值节点x0=1.00，x1=1.02，x2=1.04，x3=1.06，构造三次Lagrange插值公式，由此计算f(1.03)的近似值，并给出其实际误差；

(2) 利用插值节点x0=1，x1=1.05构造三次Hermite插值公式，由此计算f(1.03)的近似值，并给出其实际误差。

使用拉格朗日插值算法的函数lagrange.m：

%% lagrange: Lagrange插值函数  
    function Y = lagrange(X_inter,Y_inter,X)  
        % X_inter是插值节点，Y_inter是插值节点的函数值，X是待求点，Y是待求点的函数值   
        n=length(X_inter); % 插值节点个数  
        m=length(X); % 待求节点个数  
        % 下面这个循环求每个待求点的函数值  
        for X_num=1:m   
            x=X(X_num);  
            P=0.0; % P是Lagrange插值函数  
            % 下面这个循环求Lagrange插值函数P，即所有基函数L_j(x)*y_j的和  
            for j=1:n  
                L=1.0; % L是Lagrange基函数  
                % 下面这个循环求Lagrange基函数L_j(x)  
                for i=1:n  
                    if i~=j  
                        L=L*(x-X_inter(i))/(X_inter(j)-X_inter(i));  
                    end  
                end  
                P=P+L*Y_inter(j);  
            end  
            Y(X_num)=P;  
        end  
    end

使用埃尔米特插值算法的函数hermite\_interpolation.m：

%% hermite_interpolation: Hermite插值算法  
    function Y = hermite_interpolation(X_inter,Y_inter,Y_inter_de,X)  
        % Y_inter_de代表插值节点的导数值  
        n=length(X_inter); % 插值节点个数  
        m=length(X); % 待求节点个数  
        % 下面这个循环求每个待求点的函数值  
        for X_num=1:m   
            x=X(X_num);  
            H=0.0; % H是hermite插值函数  
            % 下面这个循环求hermite插值函数H  
            for j=1:n  
                L_2=1.0; % L_2是Lagrange基函数的平方  
                b=0; % b是求解H的中间量  
                % 下面这个循环求Lagrange基函数L_j(x)的平方  
                for i=1:n  
                    if i~=j  
                        L_2=L_2*((x-X_inter(i))/(X_inter(j)-X_inter(i)))^2;  
                        b=1/(X_inter(j)-X_inter(i))+b;  
                    end  
                end  
                H=H+L_2*((X_inter(j)-x)*(2*b*Y_inter(j)-Y_inter_de(j))+Y_inter(j));  
            end  
            Y(X_num)=H;  
        end  
    end

测试代码question3.m：

% 第一小题，利用Lagrange插值公式  
    syms x  
    f(x)=exp(x)*(3*x-exp(x));  
    	  
    X_inter_lag=[1.00;1.02;1.04;1.06];  
    Y_inter_lag=f(X_inter_lag);  
    X=1.03;  
    Y_lag=lagrange(X_inter_lag,Y_inter_lag,X);  
    error_lag=abs(Y_lag-f(X)); % 实际误差  
    	  
    fprintf('P(1.03)=%1.8f(cal by lagrange)\n', Y_lag);  
    fprintf('error=%1.8f(cal by lagrange)\n', error_lag);  
      
    % 第二小题，利用Hermite插值公式  
    X_inter_her=[1;1.05];  
    Y_inter_her=f(X_inter_her);  
    f_de=diff(f);  
    Y_inter_her_de=f_de(X_inter_her);  
    Y_her=hermite_interpolation(X_inter_her,Y_inter_her,Y_inter_her_de,X);  
    error_her=abs(Y_her-f(X)); % 实际误差  
    	  
    fprintf('H(1.03)=%1.8f(cal by hermite_interpolation)\n', Y_her);  
    fprintf('error=%1.8f(cal by hermite_interpolation)\n', error_her);

计算结果：

![20180522124901179][]

(25分) 利用复合求积算法计算积分，取步长为10－4

![70 6][]

使用复合梯形公式的函数compound\_trapezoidal.m：

%% compound_trapezoidal: 复合梯形公式  
    function F = compound_trapezoidal(f,a,b,h)  
        % f是被积函数，a，b是积分区间，h是步长  
        N=(b-a)/h; % 节点个数  
        f_sum=0;  
        for n=1:N-1  
            x_n=a+n*h;  
            f_sum=f_sum+f(x_n);  
        end  
        F=h/2*(f(a)+2*f_sum+f(b));  
    end

使用复合辛普森公式的函数compound\_simpson.m：

%% compound_simpson: 复合辛普森公式  
    function F = compound_simpson(f,a,b,h)  
        % f是被积函数，a，b是积分区间，h是步长  
        N=(b-a)/h; % 节点个数  
        f_sum_half=0;  
        f_sum=0;  
        x_n=a; % 保证x_n_half的起始计算  
        for n=1:N-1  
            x_n_half=x_n+h/2;  
            x_n=a+n*h;  
            f_sum_half=f_sum_half+f(x_n_half);  
            f_sum=f_sum+f(x_n);  
        end  
        F=h/6*(f(a)+4*f_sum_half+2*f_sum+f(b)); % 注意课本上x_n_half的求和是从n=0开始的  
    end

测试代码question4.m：

syms x  
    g(x)=sqrt(1+cos(x)^2);  
    a=0;  
    b=48;  
    % h=10e-4;  
    h=0.1;  
      
    % 使用复合梯形公式  
    F_compound_trapezoidal=compound_trapezoidal(g,a,b,h);  
    fprintf('cal by compound_trapezoidal: %1.8f\n', F_compound_trapezoidal);  
    int_res=int(g,[0,48]);  
    fprintf('cal by int: %1.8f\n', int_res);  % 使用matlab自带的求积分函数
    	  
    % 使用复合辛普森公式  
    F_compound_simpson=compound_simpson(g,a,b,h);  
    fprintf('cal by compound_simpson: %1.8f\n', F_compound_simpson);

计算结果：

![20180522124921914][]

其中cal by int 是使用matlab自带的求积分函数int计算出的结果。

[70]: /images/20220524/22fbd5965a554d68aea4d0c6d88beffa.png
[70 1]: /images/20220524/190284c7e9da4e10893c1bf6f3b2bfba.png
[20180522124740942]: https://img-blog.csdn.net/20180522124740942
[20180522124813813]: /images/20220524/3c2efcdfbabb4ec59a80962040b0e5fb.png
[70 2]: /images/20220524/a8c1385abef14ded8973118c7c7d1f88.png
[70 3]: /images/20220524/07651d6e9cd64178b61b27ebe87fb2b2.png
[70 4]: /images/20220524/bb1bbe9751144210a3162b5520418e23.png
[20180522124836936]: /images/20220524/85f2349b204e4ce796f39bee88f26e59.png
[70 5]: /images/20220524/ca1ad042d1a846cc8582be7dd43a7983.png
[20180522124901179]: /images/20220524/b7b4727e56954edb800e5b03ca4382b2.png
[70 6]: /images/20220524/63fca793bf734755ae15e5a79cc5a729.png
[20180522124921914]: /images/20220524/221eb5ae69f846fa9e33602f920b9552.png