快捷通道

QQ登录

只需要一步，快速开始

  注册地址  找回密码

数学建模社区-数学中国

留学>>

极限测试之Matlab与Forcal真实演练

2012-7-30 09:48| 发布者: darker50| 查看: 1317| 评论: 5|原作者: forcal

摘要: 首先要说明的是本测试系列，除了比较编译效率外，其余所有运行效率的比较都是将matlab的首次运行排出在外，因matlab程序首次运行效率较低。理论上，Forcal程序任意次运行效率都是一样的。不过话又说回来，任意程序包 ...
首先要说明的是本测试系列，除了比较编译效率外，其余所有运行效率的比较都是将matlab的首次运行排出在外，因matlab程序首次运行效率较低。理论上，Forcal程序任意次运行效率都是一样的。不过话又说回来，任意程序包含的函数，有些是需要多次运行的，而有些仅运行一次，甚至一次都不运行，故matlab函数首次运行效率较低应该是一个缺点。但如果说，matlab函数首次运行会对函数进行优化，以后运行效率会显著提高，则matlab函数首次运行效率较低就成了一个优点。

＝＝＝＝＝＝＝＝＝＝＝＝＝

本次演练要用matlab和Forcal实现两个实用函数并进行测试。当然，对脚本来说，这些函数用C/C++或Fortran来实现应是脚本的最佳选择。

＝＝＝＝＝＝＝＝＝＝＝＝＝

1、求解实系数方程组的全选主元高斯消去法：包含大量数组元素存取操作

C/C++代码：
  1. #include "stdafx.h"
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include "time.h"
  5. #include "math.h"

  7. int agaus(double *a,double *b,int n)
  8. {
  9.         int *js,l,k,i,j,is,p,q;
  10.     double d,t;
  11.     js=new int[n];
  12.     l=1;
  13.     for (k=0;k<=n-2;k++)
  14.     {
  15.                 d=0.0;
  16.         for (i=k;i<=n-1;i++)
  17.                 {
  18.           for (j=k;j<=n-1;j++)
  19.           {
  20.                           t=fabs(a[i*n+j]);
  21.               if (t>d) { d=t; js[k]=j; is=i;}
  22.           }
  23.                 }
  24.         if (d+1.0==1.0)
  25.                 {
  26.                         l=0;
  27.                 }
  28.         else
  29.         {
  30.                         if (js[k]!=k)
  31.                         {
  32.               for (i=0;i<=n-1;i++)
  33.               {
  34.                                   p=i*n+k; q=i*n+js[k];
  35.                   t=a[p]; a[p]=a[q]; a[q]=t;
  36.               }
  37.                         }
  38.             if (is!=k)
  39.             {
  40.                                 for (j=k;j<=n-1;j++)
  41.                 {
  42.                                         p=k*n+j; q=is*n+j;
  43.                     t=a[p]; a[p]=a[q]; a[q]=t;
  44.                 }
  45.                 t=b[k]; b[k]=b[is]; b[is]=t;
  46.             }
  47.         }
  48.         if (l==0)
  49.         {
  50.                         delete[] js; printf("fail\n");
  51.             return(0);
  52.         }
  53.         d=a[k*n+k];
  54.         for (j=k+1;j<=n-1;j++)
  55.         {
  56.                         p=k*n+j; a[p]=a[p]/d;
  57.                 }
  58.         b[k]=b[k]/d;
  59.         for (i=k+1;i<=n-1;i++)
  60.         {
  61.                         for (j=k+1;j<=n-1;j++)
  62.             {
  63.                                 p=i*n+j;
  64.                 a[p]=a[p]-a[i*n+k]*a[k*n+j];
  65.             }
  66.             b[i]=b[i]-a[i*n+k]*b[k];
  67.         }
  68.     }
  69.     d=a[(n-1)*n+n-1];
  70.     if (fabs(d)+1.0==1.0)
  71.     {
  72.                 delete[] js; printf("fail\n");
  73.         return(0);
  74.     }
  75.     b[n-1]=b[n-1]/d;
  76.     for (i=n-2;i>=0;i--)
  77.     {
  78.                 t=0.0;
  79.         for (j=i+1;j<=n-1;j++)
  80.                 {
  81.           t=t+a[i*n+j]*b[j];
  82.                 }
  83.         b[i]=b[i]-t;
  84.     }
  85.     js[n-1]=n-1;
  86.     for (k=n-1;k>=0;k--)
  87.         {
  88.       if (js[k]!=k)
  89.       {
  90.                   t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;
  91.           }
  92.         }
  93.     delete[] js;
  94.     return(1);
  95. }

  97.   
  98. int main(int argc, char *argv[])
  99. {
  100.         int i,j,k;
  101.     double a[4][4]=
  102.            { {0.2368,0.2471,0.2568,1.2671},
  103.              {0.1968,0.2071,1.2168,0.2271},
  104.              {0.1581,1.1675,0.1768,0.1871},
  105.              {1.1161,0.1254,0.1397,0.1490} };
  106.     double b[4]={1.8471,1.7471,1.6471,1.5471};
  107.         double aa[4][4],bb[4];
  108.         clock_t tm;

  110.         tm=clock();
  111.         for(i=0;i<10000;i++)
  112.         {
  113.                 for(j=0;j<4;j++)
  114.                 {
  115.                         for(k=0;k<4;k++)
  116.                         {
  117.                                 aa[j][k]=a[j][k];
  118.                         }
  119.                 }
  120.                 for(j=0;j<4;j++)
  121.                 {
  122.                         bb[j]=b[j];
  123.                 }
  124.                 agaus((double *)aa,bb,4);
  125.         }
  126.         printf("循环 %d 次, 耗时 %d 毫秒。\n", i,(clock()-tm));

  128.     for (i=0;i<=3;i++)
  129.         {
  130.         printf("x(%d)=%e\n",i,bb[i]);
  131.         }
  132. }
复制代码
结果：
循环 10000 次, 耗时 31 毫秒。
x(0)=1.040577e+000
x(1)=9.870508e-001
x(2)=9.350403e-001
x(3)=8.812823e-001

－－－－－－－－－

matlab 2009a代码：
  1. %file agaus.m
  2. function c=agaus(a,b,n)
  3.     js=linspace(0,0,n);
  4.     l=1;
  5.     for k=1:n-1
  6.         d=0.0;
  7.         for i=k:n
  8.           for j=k:n
  9.             t=abs(a(i,j));
  10.             if (t>d)
  11.                d=t; js(k)=j; is=i;
  12.             end
  13.           end
  14.         end
  15.         if d+1.0==1.0
  16.           l=0;
  17.         else
  18.             if js(k)~=k
  19.               for i=1:n
  20.                   t=a(i,k); a(i,k)=a(i,js(k)); a(i,js(k))=t;
  21.               end
  22.             end
  23.             if is~=k
  24.               for j=k:n
  25.                 t=a(k,j); a(k,j)=a(is,j); a(is,j)=t;
  26.               end
  27.               t=b(k); b(k)=b(is); b(is)=t;
  28.             end
  29.         end
  30.         if l==0
  31.            printf('fail\n');
  32.            c=[];
  33.            return;
  34.         end
  35.         d=a(k,k);
  36.         for j=k+1:n
  37.            a(k,j)=a(k,j)/d;
  38.         end
  39.         b(k)=b(k)/d;
  40.         for i=k+1:n
  41.           for j=k+1:n
  42.                a(i,j)=a(i,j)-a(i,k)*a(k,j);
  43.           end
  44.           b(i)=b(i)-a(i,k)*b(k);
  45.         end
  46.     end
  47.     d=a(n,n);
  48.     if abs(d)+1.0==1.0
  49.         printf('fail\n');
  50.         c=[];
  51.         return;
  52.     end
  53.     b(n)=b(n)/d;
  54.     for i=n-1:-1:1
  55.         t=0.0;
  56.         for j=i+1:n
  57.           t=t+a(i,j)*b(j);
  58.         end
  59.         b(i)=b(i)-t;
  60.     end
  61.     js(n)=n;
  62.     for k=n:-1:1
  63.       if js(k)~=k
  64.          t=b(k); b(k)=b(js(k)); b(js(k))=t;
  65.       end
  66.     end
  67.     c=b;
  68.     return;
  69. end

  71. a=[0.2368,0.2471,0.2568,1.2671;
  72.    0.1968,0.2071,1.2168,0.2271;
  73.    0.1581,1.1675,0.1768,0.1871;
  74.    1.1161,0.1254,0.1397,0.1490] ;
  75. b=[ 1.8471,1.7471,1.6471,1.5471];

  77. tic
  78. for i=1:10000
  79.     c=agaus(a,b,4);
  80. end
  81. c
  82. toc

  84. c =

  86.     1.0406    0.9871    0.9350    0.8813

  88. Elapsed time is 0.762713 seconds.
复制代码
－－－－－－－－－－

Forcal代码：
  1. !using["math","sys"];
  2. agaus(a,b,n : js,l,k,i,j,is, d,t)=
  3. {
  4.     oo{ js=array(n)},
  5.     l=1, k=0,
  6.     while{ k<n-1,
  7.         d=0.0, i=k,
  8.         while{ i<n,
  9.           j=k, while{j<n,
  10.               t=abs(a[i,j]),
  11.               if{t>d, d=t, js[k]=j, is=i},
  12.               j++
  13.           },
  14.           i++
  15.         },
  16.         which{ d+1.0==1.0, l=0,
  17.           { if{ (js[k]!=k),
  18.                 i=0, while{i<n,
  19.                   t=a[i,k], a[i,k]=a[i,js[k]], a[i,js[k]]=t,
  20.                   i++
  21.                 }
  22.             },
  23.             if{ (is!=k),
  24.                 j=k, while{j<n,
  25.                     t=a[k,j], a[k,j]=a[is,j], a[is,j]=t,
  26.                     j++
  27.                 },
  28.                 t=b[k], b[k]=b[is], b[is]=t
  29.             }
  30.           }
  31.         },
  32.         if{ (l==0),
  33.             printff("fail\r\n"),
  34.             return(0)
  35.         },
  36.         d=a[k,k],
  37.         j=k+1, while {j<n, a[k,j]=a[k,j]/d, j++},
  38.         b[k]=b[k]/d,
  39.         i=k+1, while {i<n,
  40.             j=k+1, while{j<n,
  41.                 a[i,j]=a[i,j]-a[i,k]*a[k,j],
  42.                 j++
  43.             },
  44.             b[i]=b[i]-a[i,k]*b[k],
  45.             i++
  46.         },
  47.         k++
  48.     },
  49.     d=a[(n-1),n-1],
  50.     if{ abs(d)+1.0==1.0,
  51.         printff("fail\r\n"),
  52.         return(0)
  53.     },
  54.     b[n-1]=b[n-1]/d,
  55.     i=n-2, while{i>=0,
  56.         t=0.0,
  57.         j=i+1, while{j<n, t=t+a[i,j]*b[j], j++},
  58.         b[i]=b[i]-t,
  59.         i--
  60.     },
  61.     js[n-1]=n-1,
  62.     k=n-1, while{k>=0,
  63.       if{(js[k]!=k),  t=b[k], b[k]=b[js[k]], b[js[k]]=t},
  64.       k--
  65.     },
  66.     return(1)
  67. };

  69. main(:i,a,b,aa,bb,t0)=
  70. {
  71.   oo{a=arrayinit{2,4,4 :
  72.              0.2368,0.2471,0.2568,1.2671,
  73.              0.1968,0.2071,1.2168,0.2271,
  74.              0.1581,1.1675,0.1768,0.1871,
  75.              1.1161,0.1254,0.1397,0.1490},
  76.      b=arrayinit{1,4 : 1.8471,1.7471,1.6471,1.5471},
  77.      aa=array[4,4], bb=array[4]
  78.   },
  79.   t0=clock(),
  80.   i=0, while{i<10000, aa.=a, bb.=b, agaus(aa,bb,4), i++},
  81.   outm[bb],
  82.   [clock()-t0]/1000
  83. };
复制代码
结果：
        1.04058       0.987051        0.93504       0.881282

2.125

Forcal用函数sys::A()对数组元素进行存取：
  1. !using["math","sys"];
  2. agaus(a,b,n : js,l,k,i,j,is, d,t)=
  3. {
  4.     oo{ js=array(n)},
  5.     l=1, k=0,
  6.     while{ k<n-1,
  7.         d=0.0, i=k,
  8.         while{ i<n,
  9.           j=k, while{j<n,
  10.               t=abs(A[a,i,j]),
  11.               if{t>d, d=t, A[js,k]=j, is=i},
  12.               j++
  13.           },
  14.           i++
  15.         },
  16.         which{ d+1.0==1.0, l=0,
  17.           { if{ (A[js,k]!=k),
  18.                 i=0, while{i<n,
  19.                   t=A[a,i,k], A[a,i,k]=A[a,i,A[js,k]], A[a,i,A[js,k]]=t,
  20.                   i++
  21.                 }
  22.             },
  23.             if{ (is!=k),
  24.                 j=k, while{j<n,
  25.                     t=A[a,k,j], A[a,k,j]=A[a,is,j], A[a,is,j]=t,
  26.                     j++
  27.                 },
  28.                 t=A[b,k], A[b,k]=A[b,is], A[b,is]=t
  29.             }
  30.           }
  31.         },
  32.         if{ (l==0),
  33.             printff("fail\r\n"),
  34.             return(0)
  35.         },
  36.         d=A[a,k,k],
  37.         j=k+1, while {j<n, A[a,k,j]=A[a,k,j]/d, j++},
  38.         A[b,k]=A[b,k]/d,
  39.         i=k+1, while {i<n,
  40.             j=k+1, while{j<n,
  41.                 A[a,i,j]=A[a,i,j]-A[a,i,k]*A[a,k,j],
  42.                 j++
  43.             },
  44.             A[b,i]=A[b,i]-A[a,i,k]*A[b,k],
  45.             i++
  46.         },
  47.         k++
  48.     },
  49.     d=A[a,(n-1),n-1],
  50.     if{ abs(d)+1.0==1.0,
  51.         printff("fail\r\n"),
  52.         return(0)
  53.     },
  54.     A[b,n-1]=A[b,n-1]/d,
  55.     i=n-2, while{i>=0,
  56.         t=0.0,
  57.         j=i+1, while{j<n, t=t+A[a,i,j]*A[b,j], j++},
  58.         A[b,i]=A[b,i]-t,
  59.         i--
  60.     },
  61.     A[js,n-1]=n-1,
  62.     k=n-1, while{k>=0,
  63.       if{(A[js,k]!=k),  t=A[b,k], A[b,k]=A[b,A[js,k]], A[b,A[js,k]]=t},
  64.       k--
  65.     },
  66.     return(1)
  67. };

  69. main(:i,a,b,aa,bb,t0)=
  70. {
  71.   oo{a=arrayinit{2,4,4 :
  72.              0.2368,0.2471,0.2568,1.2671,
  73.              0.1968,0.2071,1.2168,0.2271,
  74.              0.1581,1.1675,0.1768,0.1871,
  75.              1.1161,0.1254,0.1397,0.1490},
  76.      b=arrayinit{1,4 : 1.8471,1.7471,1.6471,1.5471},
  77.      aa=array[4,4], bb=array[4]
  78.   },
  79.   t0=clock(),
  80.   i=0, while{i<10000, aa.=a, bb.=b, agaus(aa,bb,4), i++},
  81.   outm[bb],
  82.   [clock()-t0]/1000
  83. };
复制代码
结果：
        1.04058       0.987051        0.93504       0.881282

1.454

－－－－－－－－－－

可以看出C/C++、matlab、Forcal耗时之比为 1 ：25：68 （Forcal不使用函数sys::A()）。
可以看出C/C++、matlab、Forcal耗时之比为 1 ：25：47 （Forcal使用函数sys::A()）。

本例Forcal耗时较长的原因在于本例程序含有大量的数组元素存取操作。

鸡蛋
鲜花
雷人
1

路过
握手

刚表态过的朋友 (1 人)

收藏 分享
发表评论

最新评论

forcal 2011-8-4 08:44
2、变步长辛卜生二重求积法：没有数组元素操作

C/C++代码：
  1. #include "stdafx.h"
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include "time.h"
  5. #include "math.h"

  7. double simp1(double x,double eps);
  8. void fsim2s(double x,double y[]);
  9. double fsim2f(double x,double y);

  11. double fsim2(double a,double b,double eps)
  12. {
  13.     int n,j;
  14.     double h,d,s1,s2,t1,x,t2,g,s,s0,ep;

  16.     n=1; h=0.5*(b-a);
  17.     d=fabs((b-a)*1.0e-06);
  18.     s1=simp1(a,eps); s2=simp1(b,eps);
  19.     t1=h*(s1+s2);
  20.     s0=1.0e+35; ep=1.0+eps;
  21.     while (((ep>=eps)&&(fabs(h)>d))||(n<16))
  22.     {
  23.                 x=a-h; t2=0.5*t1;
  24.         for (j=1;j<=n;j++)
  25.         {
  26.                         x=x+2.0*h;
  27.             g=simp1(x,eps);
  28.             t2=t2+h*g;
  29.         }
  30.         s=(4.0*t2-t1)/3.0;
  31.         ep=fabs(s-s0)/(1.0+fabs(s));
  32.         n=n+n; s0=s; t1=t2; h=h*0.5;
  33.     }
  34.     return(s);
  35. }

  37. double simp1(double x,double eps)
  38. {
  39.     int n,i;
  40.     double y[2],h,d,t1,yy,t2,g,ep,g0;

  42.     n=1;
  43.     fsim2s(x,y);
  44.     h=0.5*(y[1]-y[0]);
  45.     d=fabs(h*2.0e-06);
  46.     t1=h*(fsim2f(x,y[0])+fsim2f(x,y[1]));
  47.     ep=1.0+eps; g0=1.0e+35;
  48.     while (((ep>=eps)&&(fabs(h)>d))||(n<16))
  49.     {
  50.                 yy=y[0]-h;
  51.         t2=0.5*t1;
  52.         for (i=1;i<=n;i++)
  53.         {
  54.                         yy=yy+2.0*h;
  55.             t2=t2+h*fsim2f(x,yy);
  56.         }
  57.         g=(4.0*t2-t1)/3.0;
  58.         ep=fabs(g-g0)/(1.0+fabs(g));
  59.         n=n+n; g0=g; t1=t2; h=0.5*h;
  60.     }
  61.     return(g);
  62. }

  64. void fsim2s(double x,double y[])
  65. {
  66.         y[0]=-sqrt(1.0-x*x);
  67.     y[1]=-y[0];
  68. }

  70. double fsim2f(double x,double y)
  71. {
  72.     return exp(x*x+y*y);
  73. }

  75. int main(int argc, char *argv[])
  76. {
  77.         int i;
  78.         double a,b,eps,s;
  79.         clock_t tm;

  81.     a=0.0; b=1.0; eps=0.0001;
  82.         tm=clock();
  83.         for(i=0;i<100;i++)
  84.         {
  85.             s=fsim2(a,b,eps);
  86.         }
  87.         printf("s=%e , 耗时 %d 毫秒。\n", s, (clock()-tm));
  88. }
复制代码
结果：
s=2.698925e+000 , 耗时 78 毫秒。

－－－－－－－

matlab代码：
  1. %file fsim2.m
  2. function s=fsim2(a,b,eps)
  3.     n=1; h=0.5*(b-a);
  4.     d=abs((b-a)*1.0e-06);
  5.     s1=simp1(a,eps); s2=simp1(b,eps);
  6.     t1=h*(s1+s2);
  7.     s0=1.0e+35; ep=1.0+eps;
  8.     while ((ep>=eps)&&(abs(h)>d))||(n<16),
  9.         x=a-h; t2=0.5*t1;
  10.         for j=1:n
  11.             x=x+2.0*h;
  12.             g=simp1(x,eps);
  13.             t2=t2+h*g;
  14.         end
  15.         s=(4.0*t2-t1)/3.0;
  16.         ep=abs(s-s0)/(1.0+abs(s));
  17.         n=n+n; s0=s; t1=t2; h=h*0.5;
  18.     end
  19. end

  21. function g=simp1(x,eps)
  22.     n=1;
  23.     [y0,y1]=f2s(x);
  24.     h=0.5*(y1-y0);
  25.     d=abs(h*2.0e-06);
  26.     t1=h*(f2f(x,y0)+f2f(x,y1));
  27.     ep=1.0+eps; g0=1.0e+35;
  28.     while (((ep>=eps)&&(abs(h)>d))||(n<16))
  29.         yy=y0-h;
  30.         t2=0.5*t1;
  31.         for i=1:n
  32.             yy=yy+2.0*h;
  33.             t2=t2+h*f2f(x,yy);
  34.         end
  35.         g=(4.0*t2-t1)/3.0;
  36.         ep=abs(g-g0)/(1.0+abs(g));
  37.         n=n+n; g0=g; t1=t2; h=0.5*h;
  38.     end
  39. end

  41. %file f2s.m
  42. function [y0,y1]=f2s(x)
  43. y0=-sqrt(1.0-x*x);
  44. y1=-y0;
  45. end

  47. %file f2f.m
  48. function c=f2f(x,y)
  49.   c=exp(x*x+y*y);
  50. end

  52. %%%%%%%%%%%%%

  54. >> tic
  55. for i=1:100
  56. a=fsim2(0,1,0.0001);
  57. end
  58. a
  59. toc

  61. a =

  63.     2.6989

  65. Elapsed time is 0.995575 seconds.
复制代码
－－－－－－－

Forcal代码：
  1. fsim2s(x,y0,y1)=
  2. {
  3.   y0=-sqrt(1.0-x*x),
  4.   y1=-y0
  5. };
  6. fsim2f(x,y)=exp(x*x+y*y);
  7. //////////////////
  8. simp1(x,eps : n,i,y0,y1,h,d,t1,yy,t2,g,ep,g0)=
  9. {
  10.     n=1,
  11.     fsim2s(x,&y0,&y1),
  12.     h=0.5*(y1-y0),
  13.     d=abs(h*2.0e-06),
  14.     t1=h*(fsim2f(x,y0)+fsim2f(x,y1)),
  15.     ep=1.0+eps, g0=1.0e+35,
  16.     while {((ep>=eps)&(abs(h)>d))|(n<16),
  17.         yy=y0-h,
  18.         t2=0.5*t1,
  19.         i=1, while{i<=n,
  20.             yy=yy+2.0*h,
  21.             t2=t2+h*fsim2f(x,yy),
  22.             i++
  23.         },
  24.         g=(4.0*t2-t1)/3.0,
  25.         ep=abs(g-g0)/(1.0+abs(g)),
  26.         n=n+n, g0=g, t1=t2, h=0.5*h
  27.     },
  28.     g
  29. };

  31. fsim2(a,b,eps : n,j,h,d,s1,s2,t1,x,t2,g,s,s0,ep)=
  32. {
  33.     n=1, h=0.5*(b-a),
  34.     d=abs((b-a)*1.0e-06),
  35.     s1=simp1(a,eps), s2=simp1(b,eps),
  36.     t1=h*(s1+s2),
  37.     s0=1.0e+35, ep=1.0+eps,
  38.     while {((ep>=eps)&(abs(h)>d))|(n<16),
  39.         x=a-h, t2=0.5*t1,
  40.         j=1, while{j<=n,
  41.             x=x+2.0*h,
  42.             g=simp1(x,eps),
  43.             t2=t2+h*g,
  44.             j++
  45.         },
  46.         s=(4.0*t2-t1)/3.0,
  47.         ep=abs(s-s0)/(1.0+abs(s)),
  48.         n=n+n, s0=s, t1=t2, h=h*0.5
  49.     },
  50.     s
  51. };

  53. //////////////////

  55. mvar:
  56. t0=sys::clock(),
  57. i=0, while{i<100, a=fsim2(0,1,0.0001), i++}, a;
  58. [sys::clock()-t0]/1000;
复制代码
结果：
2.698925000624303
0.328

－－－－－－－－－

本例C/C++、matlab、Forcal运行耗时之比为 1：12.7：4.2 。

本例matlab慢的原因应该在于有大量的函数调用。另外，matlab要求每一个函数必须存为磁盘文件真的很麻烦。

本例还说明，对C/C++和Forcal脚本来说，C/C++的一条指令，Forcal平均大约用4.2条指令进行解释。
forcal 2011-8-4 09:00
3、变步长辛卜生二重求积法，写成通用的函数：没有数组元素操作

注意函数fsim2中增加了两个参数，函数句柄fsim2s用于计算二重积分时内层的上下限，函数句柄fsim2f用于计算积分函数的值。

不再给出C/C++代码，因其效率不会发生变化。

Matlab代码：
  1. %file fsim2.m
  2. function s=fsim2(a,b,eps,fsim2s,fsim2f)
  3.     n=1; h=0.5*(b-a);
  4.     d=abs((b-a)*1.0e-06);
  5.     s1=simp1(a,eps,fsim2s,fsim2f); s2=simp1(b,eps,fsim2s,fsim2f);
  6.     t1=h*(s1+s2);
  7.     s0=1.0e+35; ep=1.0+eps;
  8.     while ((ep>=eps)&&(abs(h)>d))||(n<16),
  9.         x=a-h; t2=0.5*t1;
  10.         for j=1:n
  11.             x=x+2.0*h;
  12.             g=simp1(x,eps,fsim2s,fsim2f);
  13.             t2=t2+h*g;
  14.         end
  15.         s=(4.0*t2-t1)/3.0;
  16.         ep=abs(s-s0)/(1.0+abs(s));
  17.         n=n+n; s0=s; t1=t2; h=h*0.5;
  18.     end
  19. end

  21. function g=simp1(x,eps,fsim2s,fsim2f)
  22.     n=1;
  23.     [y0,y1]=fsim2s(x);
  24.     h=0.5*(y1-y0);
  25.     d=abs(h*2.0e-06);
  26.     t1=h*(fsim2f(x,y0)+fsim2f(x,y1));
  27.     ep=1.0+eps; g0=1.0e+35;
  28.     while (((ep>=eps)&&(abs(h)>d))||(n<16))
  29.         yy=y0-h;
  30.         t2=0.5*t1;
  31.         for i=1:n
  32.             yy=yy+2.0*h;
  33.             t2=t2+h*fsim2f(x,yy);
  34.         end
  35.         g=(4.0*t2-t1)/3.0;
  36.         ep=abs(g-g0)/(1.0+abs(g));
  37.         n=n+n; g0=g; t1=t2; h=0.5*h;
  38.     end
  39. end

  41. %file f2s.m
  42. function [y0,y1]=f2s(x)
  43. y0=-sqrt(1.0-x*x);
  44. y1=-y0;
  45. end

  47. %file f2f.m
  48. function c=f2f(x,y)
  49.   c=exp(x*x+y*y);
  50. end

  52. %%%%%%%%%%%%%%%%

  54. >> tic
  55. for i=1:100
  56. a=fsim2(0,1,0.0001,@f2s,@f2f);
  57. end
  58. a
  59. toc

  61. a =

  63.     2.6989

  65. Elapsed time is 1.267014 seconds.
复制代码
－－－－－－－－

Forcal代码：
  1. simp1(x,eps,fsim2s,fsim2f : n,i,y0,y1,h,d,t1,yy,t2,g,ep,g0)=
  2. {
  3.     n=1,
  4.     fsim2s(x,&y0,&y1),
  5.     h=0.5*(y1-y0),
  6.     d=abs(h*2.0e-06),
  7.     t1=h*(fsim2f(x,y0)+fsim2f(x,y1)),
  8.     ep=1.0+eps, g0=1.0e+35,
  9.     while {((ep>=eps)&(abs(h)>d))|(n<16),
  10.         yy=y0-h,
  11.         t2=0.5*t1,
  12.         i=1, while{i<=n,
  13.             yy=yy+2.0*h,
  14.             t2=t2+h*fsim2f(x,yy),
  15.             i++
  16.         },
  17.         g=(4.0*t2-t1)/3.0,
  18.         ep=abs(g-g0)/(1.0+abs(g)),
  19.         n=n+n, g0=g, t1=t2, h=0.5*h
  20.     },
  21.     g
  22. };

  24. fsim2(a,b,eps,fsim2s,fsim2f : n,j,h,d,s1,s2,t1,x,t2,g,s,s0,ep)=
  25. {
  26.     n=1, h=0.5*(b-a),
  27.     d=abs((b-a)*1.0e-06),
  28.     s1=simp1(a,eps,fsim2s,fsim2f), s2=simp1(b,eps,fsim2s,fsim2f),
  29.     t1=h*(s1+s2),
  30.     s0=1.0e+35, ep=1.0+eps,
  31.     while {((ep>=eps)&(abs(h)>d))|(n<16),
  32.         x=a-h, t2=0.5*t1,
  33.         j=1, while{j<=n,
  34.             x=x+2.0*h,
  35.             g=simp1(x,eps,fsim2s,fsim2f),
  36.             t2=t2+h*g,
  37.             j++
  38.         },
  39.         s=(4.0*t2-t1)/3.0,
  40.         ep=abs(s-s0)/(1.0+abs(s)),
  41.         n=n+n, s0=s, t1=t2, h=h*0.5
  42.     },
  43.     s
  44. };

  46. //////////////////

  48. f2s(x,y0,y1)=
  49. {
  50.   y0=-sqrt(1.0-x*x),
  51.   y1=-y0
  52. };
  53. f2f(x,y)=exp(x*x+y*y);

  55. mvar:
  56. t0=sys::clock(),
  57. i=0, while{i<100, a=fsim2(0,1,0.0001,HFor("f2s"),HFor("f2f")), i++}, a;
  58. [sys::clock()-t0]/1000;
复制代码
结果：
2.698925000624303
0.844

－－－－－－－－

本例matlab与Forcal耗时之比为1.267014 ：0.844。Forcal仍有优势。

本例Forcal耗时增加的原因：在函数fsim2及simp1中要动态查找函数句柄fsim2s,fsim2f，并验证其是否有效，故效率下降了。
sxjm567 2012-7-30 01:48
百年不遇的好帖子，不得不顶
1055486706 2012-9-2 16:31
确实很不错
1055486706 2012-9-5 09:29
好深奥哦，肿么办？

查看全部评论(5)

相关分类

qq
收缩

  • 电话咨询

  • 04714969085

关于我们| 联系我们| 诚征英才| 对外合作| 产品服务| QQ

手机版|Archiver| |繁體中文 手机客户端  

蒙公网安备 15010502000194号

Powered by Discuz! X2.5   © 2001-2013 数学建模网-数学中国 ( 蒙ICP备14002410号-3 蒙BBS备-0002号 )     论坛法律顾问：王兆丰

GMT+8, 2024-4-27 01:25 , Processed in 0.591741 second(s), 33 queries .

回顶部