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

$\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}$

$\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)迭代法求方程组的根

$\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}}.$

$\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}}$

#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)迭代法求方程组的根

$\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)})$

#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})$

#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. 二分法求方程的根

#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. 迭代法求方程的根

#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.牛顿迭代法

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

#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. 弦割法

#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");



## 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.左边为啥能约等于右边呢？

$\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$

### 1.2.那么复化（composite）是什么意思呢 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}$

## 2.复化辛普森求积分

#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.复化梯形求积分


#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;
}
}