互联网公司职工按 996 甚至每月休 2 天的频率来开展工作,因身体和精神原因出事频繁,尤其某个“x多多”近期频繁上新闻。 最近读了马克思的《1844年经济学哲学手稿》,认识到某地还是个“初级”的“资本主义”社会。简直是天天抓理论学习,但永远都是学一套、说一套、做一套。

在《1844年经济学哲学手稿》中马克思提出了异化劳动的观点。异化的最初含义可以理解为“自我外化为非我,从而使原来与自我同一的东西,变成异己的东西”。例如,我创作一件木雕作品,把自己的思想和手艺固化到这个作品里,也就是将自己的“部分思想”和“手艺”同自己分离,形成一件雕塑,变成不是自己的的东西。马克思的观点认为:异化作为社会现象,是与阶级一起产生的,是人的物质生产机器产品变成异己力量,反过来统治人的社会现象。私有制是异化产生的社会根源,社会分工固定是它的自然缘由。在异化活动中,人的能动性丧失了。遭到异己的物质力量或精神力量的奴役,从而使人的个性不能全面发展,只能片面发展,甚至畸形发展。在资本主义社会,异化达到最严重的程度。

资本主义的工厂中,让工人延长上班时间,生产更多的产品,产生更多的利润。但产品的利润归资本家所有,反过来工人只能通过自己劳动所的工资来购买产品。工人为了购买必需品,而进行劳动,这时候,劳动就不是工人“自由自觉”的活动,具有强制性、外在性和异己性。这同互联网行业的996,以及某些实验室的研究生培养,没有任何区别。互联网行业里,大家每天回家就睡个觉,其余时间都在通勤、搬砖、吃喝拉撒,自己的产品获得的利润归资本家所有,而自己只能通过获得的工资来购买产品,期权这些等着解禁还不知要熬过多少年,前提是身体健康情况允许你继续熬。这也是互联网行业35岁以上不要的原因,当不了电池继续发电,自然要淘汰。一些打卡严格的高等教育实验室也是严格规定学生的工作时间,本来收入就只能够吃饭,还可能时常受到精神打击。在这些异化活动里,人的能动性完全丧失了,绝大部分人只能片面的发展。纵观世界,东亚三国,两个资本主义国家、一个社会主义国家,都差不多这样。还真说不好这么“卷”是因为我们自身的发达程度不够,还是人口原因或者文化原因。

马克思根据唯物主义辩证法,唯物史观和剩余价值学说,还论证过资本主义必然灭亡,我不知他看到现今的情况,会不会发表什么新的看法,毕竟他老家的工人待遇比这边的好太多。在发达国家里比较的话,估计 Musk 的公司要排除在外,我觉得他家公司的工作时长应该和我们的 996 有一拼。


一阶常微分方程初值问题,一般表示如下所示:

\[\begin{cases} y^{'}=f(x,y)\\ y(x_{0})=y_{0} \end{cases}\]

欧拉法,改进的欧拉法,龙格-库塔法都是基于同样的原理,即用切线\(y^{'}=f(x,y)\)去逼近原方程的曲线\(y=f(x)\)1

Euler_Method

那么怎么作出切线呢,\(y^{'}=f(x,y)\)就是这个切线的方程。

欧拉法计算时,令每次前进的步长为\(h\),每次根据步长h根据切线\(f(x,y)\)的方向求下一个\(y\)的位置。

\[\begin{cases} y_{i+1}=y_{i}+hK_{1}\\ K_{1}=f(x_{i},y_{i}) \end{cases}\]

改进的欧拉法,就是在欧拉法的基础上修改方向(斜率),是得每步计算的切线更贴近原曲线。

\[\begin{cases} y_{i+1}=y_{i}+\frac{h}{2}(K_{1}+K_{2})\\ K_{1}=f(x_{i},y_{i})\\ K_{2}=f(x_{i+1},y_{i}+hK_{1}) \end{cases}\]

如果设法在\([x_{i},x_{i+1}]\)内多预报几个点的斜率值,然后将它们加权平均作为平均斜率\(K^{*}\),则有可能构造出更高精度的计算格式,龙格-库塔法就是在上述改进欧拉法的基础上,继续构造新的\(K\),来达到更加精确的逼近原曲线的目的。

一般常用的方法是四阶龙格-库塔法,计算公式如下:

\[\begin{cases} K_{1}=f(x,y)\\ K_{2}=f(x_{i+\frac{1}{2}},y_{i}+\frac{h}{2}K_{1})\\ K_{3}=f(x_{i+\frac{1}{2}},y_{i}+\frac{h}{2}K_{2})\\ K_{4}=f(x_{i+1},y_{i}+hK_{3})\\ y_{i+1}=y_{i}+\frac{h}{6}(K_{1}+2K_{2}+2K_{3}+K_{4}) \end{cases}\]

四阶龙格-库塔法程序如下所示。

#include"stdio.h"
#include"stdlib.h"
void fun1(double a,double b,double h)
{  double x,y=0,i,k1,k2,k3,k4;
   for(i=1;i<=((b-a)/h)+1;i++)
 { x=i*h;
   k1=h*(1-y);
   k2=h*(1-(y+1/2.0*k1));
   k3=h*(1-(y+1/2.0*k2));
   k4=h*(1-(y+k3));
   y=y+1/6.0*(k1+2*k2+2*k3+k4);
 }  
 printf("%lf\n",y);
      }
      
void fun2(double a,double b,double h)
{  double x,y=1,i,k1,k2,k3,k4;   
   for(i=0;i<=((b-a)/h);i++)
 { x=i*h;
   k1=h*(x*y*y);
   k2=h*((x+h/2)*(y+1/2.0*k1)*(y+1/2.0*k1));
   k3=h*((x+h/2)*(y+1/2.0*k2)*(y+1/2.0*k2));
   k4=h*((x+h)*(y+k3)*(y+k3));
   y=y+1/6.0*(k1+2*k2+2*k3+k4);
 }  
 printf("%lf\n",y);
      }

int main()
{double a=0,b=1,h=0.1;
      fun1(a,b,h);
      fun2(a,b,h);
      }

1. 雅克比(Jacobi)迭代法求方程组的根

雅克比迭代法非常简单,对于一个给定的n×n方程组\(\displaystyle A\mathbf {x} =\mathbf {b}\)

\[\displaystyle A={\begin{bmatrix}a_{11}a_{12}\cdots a_{1n}\\a_{21}a_{22}\cdots a_{2n}\\\vdots \vdots \ddots \vdots \\a_{n1}a_{n2}\cdots a_{nn}\end{bmatrix}},\qquad \mathbf {x} ={\begin{bmatrix}x_{1}\\x_{2}\\\vdots \\x_{n}\end{bmatrix}},\qquad \mathbf {b} ={\begin{bmatrix}b_{1}\\b_{2}\\\vdots \\b_{n}\end{bmatrix}}.\]

可以把 A分解成两个矩阵

\[\displaystyle A=D+R\qquad \qquad D={\begin{bmatrix}a_{11}0\cdots 0\\0a_{22}\cdots 0\\\vdots \vdots \ddots \vdots \\00\cdots a_{nn}\end{bmatrix}},\qquad R={\begin{bmatrix}0a_{12}\cdots a_{1n}\\a_{21}0\cdots a_{2n}\\\vdots \vdots \ddots \vdots \\a_{n1}a_{n2}\cdots 0\end{bmatrix}}\]

方程组可以改写为 \(\displaystyle D\mathbf {x} =\mathbf {b} -R\mathbf {x}\),该公式即为迭代公式。

又可写成 \(\displaystyle \mathbf {x} ^{(k+1)}=D^{-1}(\mathbf {b} -R\mathbf {x} ^{(k)})\)。

注意这种解方程组的问题,矩阵都要满一定的收敛条件,在这里是迭代矩阵的谱半径小于1,\(\rho (D^{-1}R)<1\)。

#include"math.h"
#include"stdio.h"
#include"stdlib.h"
void main()
{
	double x1,x2,x3,x11,x22,x33;
    int i=0;
    printf("No.0008雅克比迭代法求方程组的根\n\n输入初值:");
    scanf("%lf %lf %lf",&x1,&x2,&x3);
    while(fabs(x11-(1-x2-x3)/(-8.0))>0.001 && fabs(x22-(16-x1-x3)/(-5.0))>0.001 && fabs(x33-(7-x1-x2)/(-4.0))>0.001)
  {
	x11=(1-x2-x3)/(-8.0);//写迭代式 
	x22=(16-x1-x3)/(-5.0);//写迭代式
	x33=(7-x1-x2)/(-4.0);//写迭代式
	x1=x11;
	x2=x22;
	x3=x33;
    i++;
  }	
	printf("解为:\nx1=%lf\nx2=%lf \nx3=%lf\n迭代次数:%d\n",x1,x2,x3,i);
system("pause");
} 

2. 高斯-赛德尔(Gauss-Seidel)迭代法求方程组的根

高斯-赛德尔迭代法和雅克比迭代法的区别在于,将A分解为一个上三角矩阵和下三角矩阵的形式。

\[\displaystyle A=L_{*}+U\qquad {\text{其中}}\qquad L_{*}={\begin{bmatrix}a_{11}0\cdots 0\\a_{21}a_{22}\cdots 0\\\vdots \vdots \ddots \vdots \\a_{n1}a_{n2}\cdots a_{nn}\end{bmatrix}},\quad U={\begin{bmatrix}0a_{12}\cdots a_{1n}\\00\cdots a_{2n}\\\vdots \vdots \ddots \vdots \\00\cdots 0\end{bmatrix}}\]

这样,代入迭代式后为:

\[\displaystyle L_{*}\mathbf {x} =\mathbf {b} -U\mathbf {x}\] \[\displaystyle \mathbf {x} ^{(k+1)}=L_{*}^{-1}(\mathbf {b} -U\mathbf {x} ^{(k)})\]

系数矩阵 A 严格对角占优或对称正定时,高斯-赛德尔迭代法必收敛。

#include"math.h"
#include"stdio.h"
#include"stdlib.h"
main()
{
	double x1,x2,x3,x11,x22,x33;
    int i=0;
    printf("No.0009高斯赛德尔(Gauss-Seidel)迭代法求方程组的根\n\n输入初值:");
    scanf("%lf %lf %lf",&x1,&x2,&x3);
    while(fabs(x11-(1-x2-x3)/(-8.0))>0.01 && fabs(x22-(16-x1-x3)/(-5.0))>0.01 && fabs(x33-(7-x1-x2)/(-4.0))>0.01)
  {
	x11=x1;x22=x2;x33=x3;
	x1=(1-x2-x3)/(-8.0);//写迭代式 
	x2=(16-x1-x3)/(-5.0);//写迭代式
	x3=(7-x1-x2)/(-4.0);//写迭代式
    i++;
  }	
	printf("解为:\nx1=%lf\nx2=%lf \nx3=%lf\n迭代次数:%d\n",x1,x2,x3,i);
system("pause");
}

3. 牛顿迭代法求方程组的根

这个根牛顿迭代法求解方程的根一样。。。

\[\displaystyle x_{n+1}=x_{n}-J_{F}(x_{n})^{-1}F(x_{n})\]

其中 \(J_{F}(x)\)是雅克比矩阵。

#include"math.h"
#include"stdio.h"
#include"stdlib.h"
main()
{
	double x1,x2,x11,x22;
    int i=0;
    printf("No.0010牛顿迭代法求方程组的根\n\n输入初值:");
    scanf("%lf %lf",&x1,&x2);
	x11=x1-(2*x2/(2*x2-8*x1)*(x1+2*x2-3)-2/(2*x2-8*x1)*(2*x1*x1+x2*x2-5));//写迭代式 
	x22=x2-(-4*x1/(2*x2-8*x1)*(x1+2*x2-3)+1/(2*x2-8*x1)*(2*x1*x1+x2*x2-5));
    i++;//写迭代式
 while(fabs(x11-x1)>0.001 && fabs(x22-x2)>0.001)
  { x1=x11;
	x2=x22;
	x11=x1-(2*x2/(2*x2-8*x1)*(x1+2*x2-3)-2/(2*x2-8*x1)*(2*x1*x1+x2*x2-5));//写迭代式 
	x22=x2-(-4*x1/(2*x2-8*x1)*(x1+2*x2-3)+1/(2*x2-8*x1)*(2*x1*x1+x2*x2-5));//写迭代式
    i++;
  }	
	printf("解为:\nx1=%lf\nx2=%lf \n迭代次数:%d\n",x11,x22,i);
system("pause");
} 

1. 二分法求方程的根

二分法可行的条件是有介值定理支持。

介值定理说明在实数范围内,\([a,b]\)区间上可以画出一个连续曲线。\(f(a)<u<f(b)\) 存在,那么存在\(c\), \(a<c<b\),使得\(f(c)=u\)。

二分法就是\(a\)与\(b\)异号时,介值定理的特殊情况,即存在实数\(c\)使得\(f(c)=0\)存在。

#include"math.h"
#include"stdio.h"
#include"stdlib.h"
double fun(double k)
{return 1-k-sin(k);//所求方程 
}
void main()
{double a,b,x;int i=0;
printf("No.0004二分法求方程的根\n\n输入区间左端点:");
scanf("%lf",&a);
printf("输入区间右端点:"); 
scanf("%lf",&b);
x=(a+b)/2.0;
while(/*fabs(fun(x))>0.00005||*/fabs(a-b)>0.00005)//近似范围 
{
if(fun(a)*fun(x)<0) b=x;
if(fun(b)*fun(x)<0) a=x;
x=(a+b)/2.0;
i++;}
printf("方程的根:%lf\n",x);
printf("迭代次数:%d\n\n",i);
system("pause");
}

2. 迭代法求方程的根

迭代法主要是将\(f(x)=0\)这种等式,该写成\(x=\varphi(x)\),然后给定一个初始的\(x_0\),根据公式\(x_1=\varphi(x_0)\),求出\(x_1\),接下来,比较\(x_1\)和\(x_0\)的差值大小,如果差值等于0或者一个非常非常小的数,那么则求出最终的\(x\)即为\(x_1\);如果差值还很大,则进行迭代步骤,将\(x_1\)带入到公式右侧\(\varphi(x_1)\),求出左侧的\(x_2\),重复上述比较步骤,直到\(x\)收敛于一个定值。

#include"math.h"
#include"stdio.h"
#include"stdlib.h"
double fun(double k)
{
	return pow(1+k*k,1/3.0);//所求方程 
}
void main()
{double a,b,x,x1,l=1;
int i=0,p;
printf("No.0005迭代法方程的根\n\n选择输入数据种类:\n1.数值\n2.区间\n");
printf("\n");
while(l==1)
{scanf("%d",&p);
if(p!=1&&p!=2) printf("请输入指定序号\n\n");//判断输入的序号是否正确 
else l=0;
}
switch(p)
 {
 case 1:
		printf("数据种类\n1.数值\n输入数值:");
		scanf("%lf",&x);
		x1=fun(x);i=i+1;
		break;
 case 2:printf("数据种类\n2.区间\n");
		printf("输入区间左端点:");
		scanf("%lf",&x);
		printf("输入区间右端点:");
		scanf("%lf",&x1);x1=fun((x+x1)/2.);i=i+1;
		break;
 }
while(fabs(x-x1)>0.00005)
{
x=fun(x1);
i=i+1;
x1=fun(x);
i=i+1;
}
printf("方程的根:%lf\n",x1); 
printf("迭代次数:%d\n\n",i);
system("pause");
}

3.牛顿迭代法

牛顿迭代法又称为牛顿-拉弗森方法(Newton-Raphson method),它比一般的迭代法有更高的收敛速度。

Newton-Raphson method

\(f(x_0)\)点的切线是\(f(x)\)的线性逼近。离\(f(x_0)\)点距离越近,这种逼近的效果也就越好,也就是说,切线与曲线之间的误差越小。所以我们可以说在\(f(x_0)\)点附近,\(切线\approx f(x)\)。

牛顿-拉弗森方法提出来的思路就是利用切线是曲线的线性逼近这个思想,随便找一个曲线上的\(f(x_0)\)点(为什么随便找,根据切线是切点附近的曲线的近似,应该在根点附近找,但是很显然我们现在还不知道根点在哪里),做一个切线,切线的根\(x_1\)(就是和x轴的交点)与曲线的根,还有一定的距离。我们从这个切线的根\(x_1\)出发,做一根垂线,和曲线相交于\(f(x_1)\)点,继续重复刚才的工作,多次迭代后会越来越接近曲线的根,即迭代收敛。

在计算时,首先选择一个\(x_0\),然后计算相应的\(f(x_0)\)和切线斜率\(f^{'}(x_0)\),接下来计算穿过点(\(x_0\),\(f(x_0)\))并且斜率为\(f^{'}(x_0)\)的直线和\(x\)轴的交点的\(x_1\)坐标,公式为: \(0=(x_1-x_0)f^{'}(x_0)+f(x_0)\), 将其转化为\(x_1=\varphi(x_0)\)的形式,即有:\(x_{n+1}=x_n-\frac{f(x_n)}{f^{'}(x_n)}\)

将初始的\(x_0\)代入到上述公式中,经过多次迭代,可以求得最后的根。

#include"math.h"
#include"stdio.h"
#include"stdlib.h"
double fun(double k)
{
	return k*k*k-k*k-1;//所求方程 
}
double fun1(double k)
{  
	return 3*k*k-2*k;//所求方程的一阶导函数 
}
void main()
{double a,b,x,x1,l=1;
int i=0,p;
printf("No.0006牛顿迭代法求方程的根\n\n选择输入数据种类:\n1.数值\n2.区间\n");
printf("\n");
while(l==1)
{scanf("%d",&p);
if(p!=1&&p!=2) printf("请输入指定序号\n\n");//判断输入的序号是否正确 
else l=0;
}
switch(p)
{
 case 1:
		printf("数据种类\n1.数值\n输入数值:");
		scanf("%lf",&x1);
		break;
 case 2:printf("数据种类\n2.区间\n");
		printf("输入区间左端点:");
		scanf("%lf",&x);
		printf("输入区间右端点:");
		scanf("%lf",&x1);x1=(x+x1)/2.;
		break;
 }
while(fabs(fun(x1)/fun1(x1))>0.00005)
{
x=x1-fun(x1)/fun1(x1);
x1=x;
i++;
}
printf("方程的根:%lf\n",x1); 
printf("迭代次数:%d\n\n",i);
system("pause");
}

4. 弦割法

弦割法(Secant Method)是基于牛顿法的一种改进,基本思想是用弦的斜率近似代替目标函数的切线斜率,并用割线与横轴交点的横坐标作为方程式的根的近似。

SecantMethod

公式就不推了,参见维基百科或者百度百科。

#include"math.h"
#include"stdio.h"
#include"stdlib.h"
double fun(double k)
{
	return k*k*k-k*k-1;//所求方程 
}
main()
{double a,b,x,x1,l=1,x0,x2;
int i=0,p;
printf("No.0007弦割法求方程的根\n\n选择输入种类:\n1.单点\n2.双点\n");
printf("\n");
while(l==1)
{scanf("%d",&p);
if(p!=1&&p!=2) printf("请输入指定序号\n\n");//判断输入的序号是否正确 
else l=0;
}
switch(p)
{
 case 1:
		printf("种类\n1.单点\n");
		printf("输入点1:");
		scanf("%lf",&x);
		printf("输入点2:");
		scanf("%lf",&x1);
        while(fabs(fun(x1)*(x1-x)/(fun(x1)-fun(x)))>0.00005)
        {
        x2=x1-fun(x1)*(x1-x)/(fun(x1)-fun(x));
        x1=x2;
        i++;
        }
        printf("方程的根:%lf\n",x1); 
        printf("迭代次数:%d\n\n",i);
		break;
 case 2:printf("种类\n2.双点\n");
		printf("输入点1:");
		scanf("%lf",&x);
		printf("输入点2:");
		scanf("%lf",&x1);
        while(fabs(fun(x1)*(x1-x)/(fun(x1)-fun(x)))>0.00005)
        {
        x2=x1-fun(x1)*(x1-x)/(fun(x1)-fun(x));
        x=x1;
        x1=x2;
        i++;
        }
        printf("方程的根:%lf\n",x1); 
        printf("迭代次数:%d\n\n",i);
	    break;
}
system("pause");


开始介绍本科时期学习的计算方法的内容,即求积分,求方程的根(普通的x元x次方程,方程组),所涉及到的基本步骤都是迭代循环之类的。

这个部分的程序都是C语言(极少部份C++语法)写的,system('pause')在linux下不管用。不知道为什么那时候写程序{都会换行,现在看起来好不习惯。

1.牛顿-柯特斯公式 (Newton–Cotes formulas)

后面要介绍的复化辛普森法和复化梯形法都是牛顿-柯特斯公式的特殊形式,所以先介绍一下牛顿-柯特斯公式的形式:

\[\int_a^b f(x)\mathrm{d}x\approx\sum_{i=1}^{n-1}\omega_{i}f(x_i)\]

1.1.左边为啥能约等于右边呢?

这是根据拉格朗日插值法进行的推导计算。对于拉格朗日插值法的直观理解推荐看知乎马同学的回答

拉格朗日插值法构造了穿过已知的一组点的曲线的函数表达式。

那么对于一个定积分问题,我们可以简化成求解穿过a与b两个点的曲线与X轴的面积。

Trapezoidal_rule_illustration

即 \(f(x)\) 这条曲线,可以根据a与b,采用拉格朗日插值法近似表示为 \(P_1(x)=\frac{(x-x_1)}{x_0-x_1}f(x_0)+\frac{(x-x_0)}{(x_1-x_0)}f(x_1)\),其中\(x_0=a\),\(x_1=b\)

将\(P_1(x)\)带入求解定积分(这部分的具体推导可以翻墙看The Math Guy的视频)。

\[\int_a^b f(x)\mathrm{d}x\approx\frac{1}{(x_1-x_0)}\int_{x_0}^{x_1}[(x-x_0)f(x_1)-(x-x_1)f(x_0)]\mathrm{d}x\]

最后结果等于 \(\frac{x_1-x_0}{2}[f(x_1)+f(x_0)]\),将\(\frac{x_1-x_0}{2}\) 看作 \(\omega_{i}\),假设我们已知a,b以及a,b之间一系列的点,最终可以根据拉格朗日插值法得到牛顿-柯特斯公式。

上面简化成a与b两点的例子的结果就是梯形公式\(\int_a^b f(x)\mathrm{d}x\approx\frac{b-a}{2}[f(b)+f(a)]\) 简化成a,\(\frac{a+b}{2}\),b的三点的结果,就是辛普森公式 \(\int_a^b f(x)\mathrm{d}x\approx\frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]\)

1.2.那么复化(composite)是什么意思呢 1

应用高阶牛顿-科特斯公式计算积分时,会出现数值不稳定的情况([龙格现象(Runge’s phenomenon)](https://en.wikipedia.org/wiki/Runge’s_phenomenon)),而低阶公式往往因为积分步长过大使得离散误差变大,因此,为了提高求积公式的精度,可以把积分区间分成若干个子区间,在每个子区间上使用低阶求积公式,然后将结果加起来,这种方法称为复化求积法

将区间[a,b]划分为n等分,步长为\(h=\frac{b-a}{n}\),节点为\(x_i=a+ih\), \(i=1,2,...,n+1\),在每个子区间\([x_i,x_{i+1}]\)使用梯形公式得:

\(\int_a^b f(x)\mathrm{d}x \approx \sum_{i=1}^{n}\int_{x_i}^{x_{i+1}}f(x)\mathrm{d}x = \frac{h}{2}[f(a)+2\sum_{i=1}^{n-1}f(a+ih)+f(b)]\), 其中 \(h=\frac{b-a}{n}\)

根据复化梯形公式的推导,同理可得复化辛普森公式为:

\(\int_a^b f(x)\mathrm{d}x \approx \frac{h}{6}\sum_{i=1}^{n-1}[f(x_i)+4f(x_{i+\frac{1}{2}})+f(x_{i+1})]\), 其中 \(h=\frac{b-a}{n}\)

上面这个公式对于n的取值有一些条件,辛普森法则是根据三个点的位置来推定曲线的函数表达形式,这时需要要求整个区间被分割成偶数份2,即n是偶数,公式可以写成\(\int_a^b f(x)\mathrm{d}x \approx \frac{h}{6}\sum_{j=1}^{\frac{n}{2}}[f(x_{2j-2})+4f(x_{2j-1})+f(x_{2j})]\)。

2.复化辛普森求积分

终于到了程序的部分,求\(\frac{\sin(x)}{x}\) 从a到b的积分。

#include<math.h>
#include<stdio.h>
#include<stdlib.h>
void main()
{   int n;//定义节点数 
    double a,b;//定义左右节点 
    printf("No.0001复化辛普森求积分(Simpson)\n\n输入节点数:");
    scanf("%d",&n);
    n=n-1; //计算几等分数 
    printf("输入左端点:"); 
    scanf("%lf",&a); 
    printf("输入右端点:");
    scanf("%lf",&b); 
   	double pi=3.1415927;
	n=n-1; //计算几等分数 
	double	h=(b-a)/(2*n);//计算步长
	int i,j;
	double x[2*n],y[2*n],t1=0,t2=0,t;
	x[0]=a; 
    x[2*n]=b;
	y[0]=1;//需要手动修改在左端点的值 
	y[2*n]=sin(pi/2)/(pi/2);//需要手动修改在右端点的值
	for(i=1;i<=2*n-1;i++) 
    {
        if(i*h<=b)x[i]=x[0]+i*h;
        }
	for(j=1;j<=2*n-1;j++)
	{
		y[j]=sin(x[j])/x[j];//需要手动修改
     	}
	for(i=1;i<=n;i++)
	{
        t1=t1+4*y[2*i-1];
        }
	for(i=1;i<=n-1;i++)
	{
	t2=t2+2*y[2*i];
    	}
	t=h/3*(y[0]+y[2*n]+t1+t2);//求积分 
	printf("\n步长h:%.7lf\n\n积分值t:%.7lf\n\n",h,t);
	printf("输出相对应的xy值:\n");
    for(i=0;i<=2*n;i++)
	printf("y=%.7lf---x=%.7lf\n",y[i],x[i]);
   	system("pause");
}

3.复化梯形求积分

接上,仍就求\(\frac{\sin(x)}{x}\) 从a到b的积分。


#include<math.h>
#include<stdio.h>
#include<stdlib.h>
void main()
{   int n;//定义节点数 
    double a,b;//定义左右节点 
    printf("No.0002复化梯形求积分\n\n输入节点数:");
    scanf("%d",&n);
    n=n-1; //计算几等分数 
    printf("输入左端点:"); 
    scanf("%lf",&a); 
    printf("输入右端点:");
    scanf("%lf",&b); 
   	double pi=3.1415927;
	
	double	h=(b-a)/n;//计算步长
	int i,j;
	double x[n+1],y[n+1],t1=0,t;
	x[0]=a; 
    x[n]=b;
	y[0]=1;//需要手动修改在左端点的值 
	y[n]=sin(pi/2)/(pi/2);//需要手动修改在右端点的值
	for(i=1;i<=n-1;i++) 
    {
        if(i*h<=b)x[i]=x[0]+i*h;
        }
	for(j=1;j<=n-1;j++)
	{
		y[j]=sin(x[j])/x[j];//需要手动修改
     	}
	for(i=1;i<=n-1;i++)
	{
        t1=t1+y[i];
        }
	t=h*(0.5*y[0]+0.5*y[n]+t1);//求积分 
	printf("\n步长h:%.7lf\n\n积分值t:%.7lf\n\n",h,t);
	printf("输出相对应的xy值:\n");
    for(i=0;i<=n;i++)
	printf("y=%.7lf---x=%.7lf\n",y[i],x[i]);
   	system("pause");
}

其实上面的程序写得都很烂,没有参考文献中的这个好。


#include<stdio.h>
#include<math.h>
double Function(double x)//所要计算积分的函数f(x)
{
    if(x==0)//sin(x)/x在0处的取值为1
        return 1;
    else
        return sin(x)/x;
}
//复化梯形公式
double Trapz(double a,double b,int n)
{
    double h=(b-a)/n;
    double T=0;
    for(int i=1;i<n;i++)
    {
        T=T+Function(a+i*h);
    }
    T*=2;
    T=(Function(a)+Function(b)+T)*h/2;
    return T;
}
//复化辛普森公式
double MulripleSimpson(double a,double b,int n)
{
    double h=(b-a)/n;
    double T=0;
    for(int i=0;i<n;i++)
    {
        T=T+Function(a+i*h)+4*Function(a+(i+0.5)*h)+Function(a+(i+1)*h);
    }
    T=T*h/6;
    return T;
}
void main()
{
    printf("使用复化梯形公式可得:%f\n",Trapz(0,1,8));
    printf("使用复化辛普森公式可得:%f\n",MulripleSimpson(0,1,4));
}

当年为了做这个作业,我的还搞了个整合代码(网上找了的程序,修改了一下),还包含了中矩形公式。

//数值分析--数值积分公式 
#include"iostream.h"
#include"math.h"
double c[10][10];
double f(double x)
{
  double sum=0;
  if(x==0) return 1;
  sum=sin(x)/x;//计算公式 
  return sum;
}
void initcotes(double c[][10])
{
  c[1][0]=c[1][1]=0.5;
  c[2][0]=c[2][2]=1.0/6.0;c[2][1]=2.0/3.0;
  c[3][0]=c[3][3]=1.0/8.0;c[3][1]=c[3][2]=3.0/8.0;
  c[4][0]=c[4][4]=7.0/90.0;c[4][1]=c[4][3]=16.0/45.0;c[4][2]=2.0/15.0;
  c[5][0]=c[4][5]=19.0/288.0;c[5][1]=c[5][4]=25.0/96.0;c[5][2]=c[5][3]=25.0/144.0;
}
int Trapezoid(double a,double b)
{
  cout<<"梯形公式的结果:"<<(b-a)*(f(a)+f(b))/2<<endl;
  return 1;
}
int MidRect(double a,double b)
{
  cout<<"中矩形公式的结果:"<<(b-a)*f((b+a)/2)<<endl;
  return 1;
}
int NewtonCotes(double a,double b)
{ 
  int n,k;double h;
  cout<<"请输入n的值(n值最多到5):";
  cin>>n;
  if(n>=6)
  cout<<"注意n值最多到5"<<endl;
  else{ 
  h=(b-a)/double(n);
  double sum=0;
  for(k=0;k<=n;k++)
   sum+=c[n][k]*f(a+k*h);
  cout<<"牛顿-柯特斯公式的结果:"<<(b-a)*sum<<endl;}
  return 1;
}
int STrapezoid(double a,double b)
{ 
  int n,k,q;double h;
  cout<<"1--复化梯形公式"<<endl;
  cout<<"2--复化辛普森求积公式"<<endl;
  cout<<"输入你想进行的操作:";
  cin>>q;  
  cout<<"请输入n的值:";
  cin>>n;
  h=(b-a)/double(n);
  double sum=0;
  sum+=(f(a)+f(b));
  for(k=1;k<=n-1;k++)  sum+=2*f(a+k*h);
 if(q==1)
 {
  cout<<"复化梯形公式的结果:"<<(h/2)*sum<<endl;
  return 1;
 }
  for(k=0;k<n;k++)
   sum+=4*f(a+(k+0.5)*h);
  cout<<"复化辛普森求积公式的结果:"<<(h/6)*sum<<endl;
  return 1;
}
void main()
{ cout<<"No.0003数值分析--数值积分公式"<<endl;
  double a,b;
  int p;
  cout<<"请输入积分的下限:";
  cin>>a;
  cout<<"请输入积分的上限:";
  cin>>b;
  initcotes(c);
  while(1)
  {
    cout<<"0--退出"<<endl;
 cout<<"1--梯形公式"<<endl;
 cout<<"2--中矩形公式"<<endl; 
 cout<<"3--牛顿柯特斯公式:"<<endl; 
 cout<<"4--复化公式"<<endl;
 cout<<"输入你想进行的操作:";
 cin>>p;
 switch(p)
 {
 case 1:Trapezoid(a,b);break;
 case 2:MidRect(a,b);break;
 case 3:NewtonCotes(a,b);break;
 case 4:STrapezoid(a,b);break;
 }
 if(p==0)  break;
  }
}

参考资料