数值分析III--求方程组的根
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");
}
数值分析II--求方程的根
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)，它比一般的迭代法有更高的收敛速度。
\(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)是基于牛顿法的一种改进，基本思想是用弦的斜率近似代替目标函数的切线斜率，并用割线与横轴交点的横坐标作为方程式的根的近似。
公式就不推了，参见维基百科或者百度百科。
#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");
计算方法系列I - 求积分
开始介绍本科时期学习的计算方法的内容，即求积分，求方程的根（普通的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轴的面积。
即 \(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;
}
}
参考资料
计算机图形学II--基础填充几何变换
昨天没写完，今天补上后半部分。现在回想起来计算机图形学是我本科时期上的最有意思的一门课程，其他解方程如果没有联系到实际问题，实在是太枯燥了。为啥我们的本科数学教科书不能改改，从更加应用的方向讲起呢。
扫描线算法
扫描线算法(Scanline rendering, Scanline alghorithm)主要用途是填充在屏幕上显示的几何图形。这个方法就是一个点一个点、一条线一条线，像扫描一样，把一个多边形的内部填满。 要想填充多边形内部的所有像素，需要找到一种合适的规则，能够沿着一个方向，一个像素不漏地把多边形内部填满，同时不污染多边形外部。于是上世纪六十年代，人们发明了一条水平方向的扫描线，它从\(y=0\)开始，判断与多边形的交点，这些交点把扫描线分成了若干段，之后判断哪些“段”在多边形内部，哪些“段”在多边形外部，然后把内部的部分着色，完成后，令\(y=y+1\)，即扫描线上移一格，重复之前的操作，直到扫描线不再与多边形的任何部分相交。
我的这个程序里用Bresenham’s line 的方法画多边形的边，然后用扫描线算法判断哪些像素是在多边形内部。
function scanline(x,y)
%测试数据：
% x=[10 50 30]./2;y=[30 20 70]./2;
% x=[10 30 50 20]./2;y=[20 10 50 70]./2;
% x=[20 50 110 110 50 20]./5;y=[20 10 30 80 50 70]./5;
% x=[20 25 210 110 80 20 50]./5;y=[20 5 60 80 50 70 35]./5;
% x=[20 25 100 210 110 80 20 50]./5;y=[20 5 40 30 80 50 70 35]./5;
n=length(x);
kk=1;
A=[0,0];
x=[x,x(1)];
y=[y,y(1)];
for i=1:n
[a,k]=Bresenhamline(x(i),y(i),x(i+1),y(i+1));%画边
kk=kk+k;
A=[A;a];
end
A=A(2:kk,:);
m=kk-1;
y0=min(A(:,2));
y1=max(A(:,2));
yy=y0;
datayy=[inf inf];
while yy<y1
k=0;
for i=1:m
if A(i,2)==yy
k=k+1;
D(yy,k)=A(i,1);
end
end
d0=min(D(yy,1:k));
d1=max(D(yy,1:k));
for j=d0:d1-1
% pause(0.001);
% plot(j,yy,'ro');
datayy=[datayy;j yy];
end
yy=yy+1;
end
x0=min(A(:,1));
x1=max(A(:,1));
xx=x0;
dataxx=[inf inf];
while xx<x1
k=0;
for i=1:m
if A(i,1)==xx
k=k+1;
D(xx,k)=A(i,2);
end
end
d0=min(D(xx,1:k));
d1=max(D(xx,1:k));
for j=d0:d1-1
% pause(0.001);
% plot(xx,j,'ro');
dataxx=[dataxx;xx j];
end
xx=xx+1;
end
if size(dataxx(:,1))>size(dataxx(:,1))
for i=2:size(dataxx(:,1))
for j=2:size(datayy(:,1))
if dataxx(i,1)==datayy(j,1) && dataxx(i,2)==datayy(j,2)
plot(dataxx(i,1),dataxx(i,2),'ro');
pause(0.001);
end
end
end
else
for i=2:size(datayy(:,1))
for j=2:size(dataxx(:,1))
if datayy(i,1)==dataxx(j,1) && datayy(i,2)==dataxx(j,2)
plot(datayy(i,1),datayy(i,2),'ro');
pause(0.001);
end
end
end
end
end
其他图形变换
最终我为了展示自己所有的画图方法，搞了个大demo程序，把所有画线和几何图形的变换都囊括在一张图里。程序里可能还有些bug。 这其中包括：旋转变化，平移变换，比例变换，对称变换（关于x轴），错切变换，相对(2,2)点的旋转变换。
错切变换（transvection）是啥？就是把矩形变成平行四边形的变换。
function demo
figure
subplot(4,2,1)
[Dx Dy]=DDALine(4,6,8,10)%DDALine(x(1),y(1),x(2),y(2))
[X2]=Bresenhamline(1,2,6,7)%Bresenhamline(x(1),y(1),x(2),y(2))
Bx=X2(:,1);
By=X2(:,2);
%----------------------------------------------
%二维几何变换
s=45;
T=[cos(s) sin(s) 0;
-sin(s) cos(s) 0;
0 0 1];
title('旋转变换(逆时针)');
xlabel('x');
ylabel('y');
for i=1:size(Dx(:))
z=[double(Dx(i)) double(Dy(i)) 1];
z=z*T;
XX(i)=z(1);
YY(i)=z(2);
plot(XX,YY,'*- k')
end
for i=1:size(Bx(:))
z=[double(Bx(i)) double(By(i)) 1];
z=z*T;
XXX(i)=z(1);
YYX(i)=z(2);
plot(XXX,YYX,'*- y')
end
%----------------------------------------------
subplot(4,2,2)
T=[1 0 0;
0 1 0;
4 5 1];
for i=1:size(Bx(:))
z=[double(Bx(i)) double(By(i)) 1];
z=z*T;
XXX(i)=z(1);
YYX(i)=z(2);
plot(XXX,YYX,'*- g')
title('平移变换');
xlabel('x');
ylabel('y');
end
hold on
plot(Bx,By,'* b');
%----------------------------------------------
subplot(4,2,3)
T=[2 0 0;
0 2 0;
0 0 1];
for i=1:size(Bx(:))
z=[double(Bx(i)) double(By(i)) 1];
z=z*T;
XXX(i)=z(1);
YYX(i)=z(2);
plot(XXX,YYX,'*- g')
title('比例变换');
xlabel('x');
ylabel('y');
end
hold on
plot(Bx,By,'* b');
%----------------------------------------------
subplot(4,2,4)
T=[1 0 0;
0 -1 0;
0 0 1];
for i=1:size(Bx(:))
z=[double(Bx(i)) double(By(i)) 1];
z=z*T;
XXX(i)=z(1);
YYX(i)=z(2);
plot(XXX,YYX,'*- g')
title('对称变换(关于x轴)');
xlabel('x');
ylabel('y');
end
hold on
plot(Bx,By,'* b');
%----------------------------------------------
subplot(4,2,5)
T=[0 -1 0;
-1 0 0;
0 0 1];
for i=1:size(Bx(:))
z=[double(Bx(i)) double(By(i)) 1];
z=z*T;
XXX(i)=z(1);
YYX(i)=z(2);
plot(XXX,YYX,'*- g')
title('对称变换(关于y=-x轴)');
xlabel('x');
ylabel('y');
end
hold on
plot(Bx,By,'* b');
%----------------------------------------------
subplot(4,2,6)
T=[1 -1 0;
2 1 0;
0 0 1];
zx=[1 5 3 1];
zy=[3 4 0 3];
clear XXX
clear YYX
for i=1:size(zx(:))
z=[double(zx(i)) double(zy(i)) 1];
z=z*T;
XXX(i)=z(1);
YYX(i)=z(2);
plot(XXX,YYX,'*- g')
title('错切变换');
xlabel('x');
ylabel('y');
end
hold on
plot(zx,zy,'*- b');
%----------------------------------------------
subplot(4,2,7)
T1=[1 0 0;
0 1 0;
-2 -2 1];
T=[cos(s) sin(s) 0;
-sin(s) cos(s) 0;
0 0 1];
clear XXX
clear YYX
for i=1:size(Bx(:))
z=[double(Bx(i)) double(By(i)) 1];
z=z*T1*T*(-1*T1);
XXX(i)=z(1);
YYX(i)=z(2);
plot(XXX,YYX,'*- g')
title('相对(2,2)点的旋转变换');
xlabel('x');
ylabel('y');
end
hold on
plot(Bx,By,'* b');
end
树根和回文串
还是那个code文件夹里的东西，现在写代码的习惯里必须再增加一条，写明这个程序的作用。我整理了不少图形图像和online judge的程序，每次回想程序的用途都是一件非常麻烦的事情，有的时候根本就想不起来。
树根
这个问题非常的经典，具体的介绍可以参照知乎中leetcode.wang的介绍
对于给定的数字，将其每一位上的数字相加得到新的数字，一直重复这个过程，直到这个数小于10，将这个数输出。
原数是n，树根就可以表示成(n-1) mod 9 + 1。 主要的用途是计算模运算的同余，对于非常大的数字的情况下可以节省很多时间。数字根可作为一种检验计算正确性的方法。例如，两数字的和的数根等于两数字分别的数根的和。
当时做这道题的时候，想到了树根的问题，但是没有求出其表达式，因而采用了一个比较苯的办法来做计算，即排除数字9,在一个多位数字中，排除9和0的其他数字之和不影响最终的树根结果。那么我就在计算的时候默认排除了9和0的加法运算。
写了一大长串程序，其中的char - '0' 是ASCII运算，将字符转换成其对应的数字，这个程序也是本篇文章中唯一可以正常运行的程序。
#include <stdio.h>
int main()
{
unsigned long a=0,b=0;
int flag=2;
char c;
while(1)
{
a=0;b=0;flag=2;
c=getchar();
if(c=='0') break;
else
{ a=a+(c-'0');
while((c=getchar())!='\n')
{
if(((c-'0')!=9)&&((c-'0')!=0))
{
a=a+(c-'0');
}
}
// c='0';
while(flag>1)
{ if(a>9)
{
while(a!=0)
{ b=b+a%10;
a=a/10;
}
if(b>9)
{
flag=2;
a=b;
b=0;
}
else{
flag=0;
a=b;
b=0;
break;
}
}
else
{
break;
}
}
printf("%d\n",a);
}
}
}
谁知道其实答案如此简单(⊙ˍ⊙)。
public int addDigits(int num) {
return (num - 1) % 9 + 1;
}
回文串
这程序里有bug，目前无法自动停止。输入一串数字，判断这窜数字是否有回文结构。这个题目主要靠算法思想，但是我没什么思想，只会最简单的按位比较.
例如：输入123321，返回YES；输入123231,返回NO。
#include <cstring>
#include<iostream>
using namespace std;
int isPalindrome(char str[])
{
int len=strlen(str);
int i,j;
i=0;
j=len-1;
while(i<j)
{
if(str[i++]!=str[j--]) return 0;
}
return 1;
}
int main(){
char str[10001];
int n,i=0;
cin>>n;
while(i<=n-1)
{ i=i+1;
cin>>str;
int a=0;
a=isPalindrome(str);
if(a==1) cout<<"YES"<<endl;
else cout<<"NO"<<endl;
}
return 0;
}
其他
我忘记了这个程序是做什么的，程序名称是MaxArray，估计是个找最大array相关的程序。从这个例子中深刻的体会到，写代码注释的重要性。
#include <stdio.h>
//#include <stdlib.h>
#define MAX 130
int maxSum(int b[], int n)
{
int max = b[0];
int sum = 0;
for(int i = 0; i < n; ++i)
{
if(sum < 0) sum = b[i];
else sum += b[i];
if(sum > max) max = sum;
}
return max;
}
int solve(int a[][MAX], int n)
{
int b[MAX];int i,j,k,sum,tmp;
sum = a[0][0];
for(i = 0; i < n; ++i)
{
for(k = 0; k < n; ++k) b[k] = 0;
for(j = i; j < n; ++j)
{
for(k =0; k < n; ++k) b[k] += a[j][k];
if(sum < (tmp = maxSum(b,n)) ) sum = tmp;
}
}
return sum;
}
int main()
{
int arr[MAX][MAX];
int i,j,n;
while(scanf("%d",&n)!=EOF)
{
for(i = 0; i < n; ++i)
for(j = 0; j < n; ++j)
scanf("%d",&arr[i][j]);
printf("%d",solve(arr,n));
break;
}
// system("pause");
return 0;
}
