开始介绍本科时期学习的计算方法的内容,即求积分,求方程的根(普通的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;
  }
}

参考资料


昨天没写完,今天补上后半部分。现在回想起来计算机图形学是我本科时期上的最有意思的一门课程,其他解方程如果没有联系到实际问题,实在是太枯燥了。为啥我们的本科数学教科书不能改改,从更加应用的方向讲起呢。

扫描线算法

扫描线算法(Scanline rendering, Scanline alghorithm)主要用途是填充在屏幕上显示的几何图形。这个方法就是一个点一个点、一条线一条线,像扫描一样,把一个多边形的内部填满。 要想填充多边形内部的所有像素,需要找到一种合适的规则,能够沿着一个方向,一个像素不漏地把多边形内部填满,同时不污染多边形外部。于是上世纪六十年代,人们发明了一条水平方向的扫描线,它从\(y=0\)开始,判断与多边形的交点,这些交点把扫描线分成了若干段,之后判断哪些“段”在多边形内部,哪些“段”在多边形外部,然后把内部的部分着色,完成后,令\(y=y+1\),即扫描线上移一格,重复之前的操作,直到扫描线不再与多边形的任何部分相交。

ScanLine

我的这个程序里用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)是啥?就是把矩形变成平行四边形的变换。

demo

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的介绍

DigitalRoot

对于给定的数字,将其每一位上的数字相加得到新的数字,一直重复这个过程,直到这个数小于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;
}


依旧是整理code文件夹里的东西,Fortran语言基本上被我当作化石了,除了在久远的“遗产代码”中你会看到这门语言,当今,应该没有人将它用作主要的编程语言。

我在学习Fortran的时候还是用的 g77,但是现在发现,在我的系统里只能使用gfortran来编译了。

下面是写过的唯一一个代码,其中C开头的是注释,编译方法是 gfortran helloworld.f -o helloworld

传统的 Fortran 程序只能用大写字符书写,而且每行前六个字符为特定用途所保留。第一列为字符 C 所保留,用来表征整行都是注释。第二列到第六列是为标号预留的。代码从第七列开始。下面是示例程序采用的是传统的 Fortran 格式。

C   helloworld.f
C
      PROGRAM HELLOWORLD
      WRITE(*,10)
   10 FORMAT('hello, world')
      END PROGRAM HELLOWORLD
C   http://wiki.ubuntu.org.cn/Compiling_Fortran77
c   g77 helloworld.f -o helloworld

希望这些老古董能尽早退役吧。


为什么我今天这么闲,可以更水这么多篇博客,那是因为昨天在家办公感觉比在办公室还累,所以索性不干了。(*  ̄︿ ̄)

本篇内容所有程序都是MATLAB/GNU Octave上运行的。

注意:有没有最后的end可能是上面两个编译环境的区别,但是我没有安装MATLAB,所以没法测试MATLAB现在的函数组成中是否仍旧不带最后一个end

数值微分法画线

众所周知,在屏幕上连续的点组成了线条。那么如何在计算机中模拟这些由“点”生成的线呢,这就是计算机图形学要研究的内容,其中最简单的一个问题就是画直线。

数值微分法( Digital Differential Analyzer,DDA)是一种处理此种问题的经典算法。

算法思想:

  1. 给定一条线段的起点\((x_1, y_1)\)和终点\((x_2, y_2)\),分别计算在两轴上的差值 \(\Delta y = y_{2} - y_{1}\) 和 \(\Delta x = x_{1} - x_{2}\)。
  2. 比较\(\Delta y\)和\(\Delta x\)二者谁比较大,大的那一个就作为遍历的总步数。\(steps = max(\Delta x, \Delta y)\)。我的程序中简化了这个部分,直接把x作为主序方向。
  3. 计算在x和y方向上单步的步距。 \(d_x = \frac{\Delta x}{steps}\),\(d_y = \frac{\Delta y}{steps}\)。我的程序里只计算了y方向的单步距离,x方向一直是“进1”。
  4. 我的程序中根据\(x_0\)和\(x_1\)的大小决定划线的方向,y增长的部分是根据斜率和目前点的位置来计算的。这个过程涉及大量的浮点运算,效率上是比较低的。
function [Xdata Ydata]=DDALine(x0,y0,x1,y1)
% DDALine(8,9,2,6)
% DDALine(2,6,8,9)
% DDALine(1,9,6,5)
% DDALine(6,5,1,9)
deltaX=x1-x0;
deltaY=y1-y0;
if deltaX~=0&&deltaY~=0
    Slope=double(deltaY)/double(deltaX);
    CurrentY=double(y0);
    j=1;
    if x0<x1
        for CurrentX=int8(x0):int8(x1)
                plot(CurrentX,int8(CurrentY-0.5),'* r');
                Xdata(j)=CurrentX;
                Ydata(j)=int8(CurrentY-0.5);
                hold on;
                
                CurrentY=Slope+CurrentY;
                j=j+1;
        end
    end
     if x0>x1
        for CurrentX=int8(x0):-1:int8(x1)
            plot(CurrentX,int8(CurrentY-0.5),'* r');
            Xdata(j)=CurrentX;
            Ydata(j)=int8(CurrentY-0.5);
            hold on;
            
            CurrentY=-Slope+CurrentY;
            j=j+1;
        end
    end
    title('DDAline');
    xlabel('x');
    ylabel('y');
    hold on
end
if deltaX==0&&deltaY~=0
    plot(x1,min(y0,y1):max(y0,y1),'* r');
end
if deltaX~=0&&deltaY==0
    plot(min(x0,x1):max(x0,x1),y0,'* r');
end
grid on;
end

Bresenham’s line algorithm

Bresenham line’s algorithm 是DDA的一种改进算法。它与DDA相比有质量和效率的两点改进:

  1. 以y方向为例,Bresenham line 根据斜率的方向,决定每一个点是离\(y_1\)还是\(y_2\)近,从而给出下一个点的位置究竟是y+1还是y-1。
  2. 由于上述原因,该方法避免了大量的浮点小数计算。
function  [a,k]=Bresenhamline(x0,y0,x1,y1)
% Bresenhamline(8,9,2,6)
% Bresenhamline(2,6,8,9)
% Bresenhamline(1,9,6,5)
% Bresenhamline(6,5,1,9)
dx=x1-x0;
dy=y1-y0;

d1=abs(2*dx);
d2=abs(2*dy);
x=x0;y=y0;

plot(x,y,'*');
a=[x y];
k=1;
if(abs(dx)>=abs(dy))
    e=-abs(dx);
while abs(x)~=abs(x1)
    if(dx>=0) x=x+1;end
    if(dx<0) x=x-1;end
    e=e+d2;
    if e>0
        if(dy>=0)y=y+1;end
        if(dy<0)y=y-1;end
        e=e-d1;          
    end
    hold on;
    a=[a;x y];k=k+1;
    hold on;
    plot(x,y,'*');
end
end
if(abs(dx)<abs(dy))
    e=-abs(dy);
while abs(y)~=abs(y1)
    if(dy>=0) y=y+1;end
    if(dy<0) y=y-1;end
    e=e+d1;
    if e>0
        if(dx>=0) x=x+1;end
            if(dx<0) x=x-1;end
        e=e-d2;          
    end
    hold on;
    a=[a;x y];k=k+1;
    hold on;
    plot(x,y,'*');
end
end
hold on;
grid on;
end