《计算方法》上机实验报告（华科软院）

叁歲伎倆 2023-06-17 14:58 22阅读 0赞

1. (25分)计算积分

![format_png][]， n=0,1,2,…,20

若用下列两种算法

(A) ![format_png 1][]

(B) ![format_png 2][]

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

代码实现(C语言)：

#include <stdio.h>  
    #include <math.h>  
      
    void main(){  
      
        double i=0.18232155679395;//ln6-ln5  
        printf("算法A：\n");  
        printf("I(0)=0.18232155679395\n");  
        for(int n=1;n<=20;n++)  
        {  
            printf("I(%d)=%lf\n",n,i);  
            i=-5*i+pow(n,-1);  
        }  
      
        double j=0.007997523028232166;  
        printf("\n算法B：\n");  
        for(int n=20;n>=0;n--)  
        {  
            printf("I(%d)=%lf\n",n,j);  
            j=0.2*(pow(n,-1)-j);  
        }  
          
    }

运行结果：

![format_png 3][]

分析：算法A每向前推进一步，其计算值的舍入误差便增长5倍；而算法B每向后推进一步，其舍入误差便减少5倍。因此算法B更合理

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

![format_png 4][]， n≥1，x0给定

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

(2) 利用牛顿迭代求解方程

![format_png 6][]

设预定精度ε=10-10。

代码实现(MATLAB):

newton.m：

% newton: Newton迭代法  
    function [x1,k] = newton(f,f_de,varible,x0,eplison)  
        dx = -subs(f_de,varible,x0)\subs(f,varible,x0);  
        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

上机作业2.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')];  
    eplison=10e-10;  
    variable=[x;y];  
    X0=[1;1];  
      
    [x,k]=newton(F,F_de,variable,X0,eplison);  
    fprintf('root:\n', k);  
    for i=1:length(x)  
        fprintf('%1.6f\n',x(i));  
    end

运行结果：

>> 上机作业2  
    root:  
    0.724712  
    1.025858

3. (25分) 已知![format_png 7][]，

(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)的近似值，并给出其实际误差。

(1)代码实现(C语言)：

#include <stdio.h>  
    #include <math.h>  
      
    double f(double x)  
    {  
        double result =((x-1.02)*(x-1.04)*(x-1.06))/((1.00-1.02)  
                       *(1.00-1.04)*(1.00-1.06))*0.765789386+  
                       ((x-1.00)*(x-1.04)*(x-1.06))/((1.02-1.00)  
                               *(1.02-1.04)*(1.02-1.06))*0.795366778+  
                       ((x-1.00)*(x-1.02)*(x-1.06))/((1.04-1.00)  
                               *(1.04-1.02)*(1.04-1.06))*0.82268817+  
                       ((x-1.00)*(x-1.02)*(x-1.04))/((1.06-1.00)  
                               *(1.06-1.02)*(1.06-1.04))*0.84752258;  
        return result;  
    }  
      
    void main()  
    {  
        double x;  
        printf("please input x: ");  
        scanf("%lf",&x);  
        printf("f(%lf)=%lf\n",x,f(x));  
    }

运行结果：

![format_png 8][]

实际误差：0.8093236189017053-0.809324=-3.8109829469945566e-7

(2) 代码实现(C语言)：

#include <stdio.h>  
    #include <math.h>  
      
    double L1(double x){  
        return (x-1.02)*(x-1.05)/((1.00-1.02)*(1-1.05));  
    }  
    double L2(double x){  
        return (x-1.00)*(x-1.05)/((1.02-1.00)*(1.02-1.05));  
    }  
    double L3(double x){  
        return (x-1.00)*(x-1.02)/((1.05-1.00)*(1.05-1.02));  
    }  
    double a1(double x){  
        return (1-2.0*(x-1.00)*(-70.0))*pow(L1(x),2.0);  
    }  
    double a2(double x){  
        return (1-2.0*(x-1.02)*(50.0/3.0))*pow(L2(x),2.0);  
    }  
    double a3(double x){  
        return (1-2.0*(x-1.05)*(160/3))*pow(L3(x),2.0);  
    }  
    double belta1(double x){  
        return (x-1.00)*pow(L1(x),2.0);  
    }  
    double belta2(double x){  
        return (x-1.02)*pow(L2(x),2.0);  
    }  
    double belta3(double x){  
        return (x-1.05)*pow(L3(x),2.0);  
    }  
    double H3(double x){  
        double y1,y2,y3,derivative_y1,derivative_y2,derivative_y3;  
        y1=0.7657893864464859;  
        y2=0.7953667788517541;  
        y3=0.8354311093313167;  
        derivative_y1=1.5315787728929717;  
        derivative_y2=1.4243418718656515;  
        derivative_y3=1.242214550953158;  
        double result = y1*a1(x)+y2*a2(x)+y3*a3(x)+derivative_y1*belta1(x)+derivative_y2*belta2(x)+derivative_y3*belta3(x);  
        return result;  
    }  
      
    void main(){  
        double x;  
        printf("please input x: ");  
        scanf("%lf",&x);  
        printf("f(%lf)=%lf\n",x,H3(x));  
    }

运行结果：

![format_png 9][]

实际误差：0.8093236189017053-0.808878=0.000445618901705358

4. (25分) 利用复合求积算法计算积分

![format_png 10][]

取步长为10－4

代码实现(MATLAB)：

algorithm\_1.m：

%% algorithm_1: 复合梯形公式  
    function F = algorithm_1(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

algorithm\_2.m：

%% algorithm_2: 复合辛普森公式  
    function F = algorithm_2(f,a,b,h)  
        N=(b-a)/h;  
        f_sum_half=0;  
        f_sum=0;  
        x_n=a;  
        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));  
    end

上机作业4.m：

syms x  
    g(x)=sqrt(1+cos(x)^2);  
    a=0;  
    b=48;  
    % h=10e-4;  
    h=0.1;  
    F_algorithm_1=algorithm_1(g,a,b,h);  
    fprintf('algorithm_1: %1.6f\n', F_algorithm_1);  
    F_algorithm_2=algorithm_2(g,a,b,h);  
    fprintf('algorithm_2: %1.6f\n', F_algorithm_2);

运行结果：

>> 上机作业4  
    algorithm_1: 58. 470466  
    algorithm_2: 58. 462540

[format_png]: https://imgconvert.csdnimg.cn/aHR0cDovL2RvY2ltZy5uY21lbS5jb20vNDAwL2NlZWRlMDMyLnBuZw?x-oss-process=image/format,png
[format_png 1]: https://imgconvert.csdnimg.cn/aHR0cDovL2RvY2ltZy5uY21lbS5jb20vNDAwLzM3ZTcxNWZhLnBuZw?x-oss-process=image/format,png
[format_png 2]: https://imgconvert.csdnimg.cn/aHR0cDovL2RvY2ltZy5uY21lbS5jb20vNDAwLzkyODNjMTI4LnBuZw?x-oss-process=image/format,png
[format_png 3]: https://imgconvert.csdnimg.cn/aHR0cDovL2RvY2ltZy5uY21lbS5jb20vNDAwLzEwMGFmY2FkLnBuZw?x-oss-process=image/format,png
[format_png 4]: https://imgconvert.csdnimg.cn/aHR0cDovL2RvY2ltZy5uY21lbS5jb20vNDAwLzkyNmFiOGU1LnBuZw?x-oss-process=image/format,png
[format_png 5]: https://imgconvert.csdnimg.cn/aHR0cDovL2RvY2ltZy5uY21lbS5jb20vNDAwL2JmMDBiMDIwLnBuZw?x-oss-process=image/format,png
[format_png 6]: https://imgconvert.csdnimg.cn/aHR0cDovL2RvY2ltZy5uY21lbS5jb20vNDAwL2I5MjlmZGM5LnBuZw?x-oss-process=image/format,png
[format_png 7]: https://imgconvert.csdnimg.cn/aHR0cDovL2RvY2ltZy5uY21lbS5jb20vNDAwLzBjYWY4M2Y2LnBuZw?x-oss-process=image/format,png
[format_png 8]: https://imgconvert.csdnimg.cn/aHR0cDovL2RvY2ltZy5uY21lbS5jb20vNDAwLzVjZDdmOWM2LnBuZw?x-oss-process=image/format,png
[format_png 9]: https://imgconvert.csdnimg.cn/aHR0cDovL2RvY2ltZy5uY21lbS5jb20vNDAwLzY3MzFjNmU2LnBuZw?x-oss-process=image/format,png
[format_png 10]: https://imgconvert.csdnimg.cn/aHR0cDovL2RvY2ltZy5uY21lbS5jb20vNDAwLzk1NjkyMjZhLnBuZw?x-oss-process=image/format,png