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

1.     (25分)计算积分

，  n=0,1,2,…,20

若用下列两种算法

(A)

(B)

试依据积分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);  
    }  
      
}  

运行结果：

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

 

 

 

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

，  n≥1，x0给定

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

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

        设预定精度ε=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分) 已知，

(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));  
}  

 

运行结果：

实际误差：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));  
}  

运行结果：

实际误差：0.8093236189017053-0.808878=0.000445618901705358

 

 

 

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

取步长为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