歡迎訪問 Forcal程序設計

FcMath庫計算舉例

例子1 矩陣輸出與子矩陣:

!using["math"];       //使用命名空間math
main(:a)=oo{         
//一般在oo函數中調用FcMath函數
    a=rand[6,5],     
//生成6×5矩陣a,用0~1之間隨機數初始化
    a.outm(),        
//輸出矩陣a
    a(neg:3).outm(), 
//取矩陣a第4列所有元素組成子矩陣,并輸出
    a(3:neg).outm(), 
//取矩陣a第4行所有元素組成子矩陣,并輸出
    a(3,5:2,3).outm()
//取矩陣a第4~6行,3~4列所有元素組成子矩陣,并輸出
};

    結果:

0.211319     4.91638e-002 0.144638     0.153259     0.852615
0.630646     0.927048     0.440308     0.162857     0.556854
0.43309      0.34552      0.563919     0.937164     0.209641
0.603271     0.727676     0.130951     5.35736e-002 0.197937
0.576004     0.747589     1.17645e-002 0.363892     0.280777
0.646454     0.381088     0.58551      0.26387      0.93692

0.153259
0.162857
0.937164
5.35736e-002
0.363892
0.26387

0.603271     0.727676     0.130951     5.35736e-002 0.197937

0.130951     5.35736e-002
1.17645e-002 0.363892
0.58551      0.26387

例子2 效率測試:

    Matlab代碼:

clear all
clc
tic
k = zeros(5,5);
% //生成5×5全0矩陣
% 循環計算以下程序段100000次:
for m = 1:100000
a = rand(5,7);
b = rand(7,5); 
%//生成5×7矩陣a,7×5矩陣b,用0~1之間的隨機數初始化
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
end
k
toc

    Forcal代碼:

!using["math","sys"];
(:t0,k,i,a,b)=
{
    t0=clock(),
    oo{k=zeros[5,5]},
    i=0,(i<100000).while{
        oo{
            a=rand[5,7], b=rand[7,5],
            k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
        },
        i++
    },
    k.outm(),
    [clock()-t0]/1000
};

    結果:

274943 274875 274944 274988 275039
274911 275007 275034 275057 274952
274915 274963 274943 274935 274972
274936 275000 275043 275027 275084
274940 274946 274999 274969 274973

3.484  秒

例子3 :ndgrid函數及多維數組的輸出

!using["math"];
mvar:
oo{
    ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],
    a1.outm[5,1,1],
    a2.outm[5,1,1],
    a3.outm[5,1,1],
    a4.outm[5,1,1],
    a=a1+a2+a3+a4,
    a.outm[5,1,1],
    Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]
};

    說明:

    linspace(8,12,5):生成一維數組,共5個元素8~12
    a1.outm[5,1,1]:輸出**數組a1,連下標一起輸出
    Sum[a]:設有m維數組(含矩陣):a(n1,n2,... ...,nm),則Sum(a,i)對第i維求和,返回一個m-1維數組。若a是一維數組、1×k矩陣、k×1矩陣、i<1或者i>m,則Sum函數返回所有數組元素的和。Sum(a)相當于Sum(a,m)。

    結果(最終求和結果是1320):

(0,0,0,*)              1.            1.            1.            1.            1.
(0,0,1,*)              1.            1.            1.            1.            1.
(0,1,0,*)              1.            1.            1.            1.            1.
(0,1,1,*)              1.            1.            1.            1.            1.
(0,2,0,*)              1.            1.            1.            1.            1.
(0,2,1,*)              1.            1.            1.            1.            1.
(1,0,0,*)              2.            2.            2.            2.            2.
(1,0,1,*)              2.            2.            2.            2.            2.
(1,1,0,*)              2.            2.            2.            2.            2.
(1,1,1,*)              2.            2.            2.            2.            2.
(1,2,0,*)              2.            2.            2.            2.            2.
(1,2,1,*)              2.            2.            2.            2.            2.

(0,0,0,*)              3.            3.            3.            3.            3.
(0,0,1,*)              3.            3.            3.            3.            3.
(0,1,0,*)              4.            4.            4.            4.            4.
(0,1,1,*)              4.            4.            4.            4.            4.
(0,2,0,*)              5.            5.            5.            5.            5.
(0,2,1,*)              5.            5.            5.            5.            5.
(1,0,0,*)              3.            3.            3.            3.            3.
(1,0,1,*)              3.            3.            3.            3.            3.
(1,1,0,*)              4.            4.            4.            4.            4.
(1,1,1,*)              4.            4.            4.            4.            4.
(1,2,0,*)              5.            5.            5.            5.            5.
(1,2,1,*)              5.            5.            5.            5.            5.

(0,0,0,*)              6.            6.            6.            6.            6.
(0,0,1,*)              7.            7.            7.            7.            7.
(0,1,0,*)              6.            6.            6.            6.            6.
(0,1,1,*)              7.            7.            7.            7.            7.
(0,2,0,*)              6.            6.            6.            6.            6.
(0,2,1,*)              7.            7.            7.            7.            7.
(1,0,0,*)              6.            6.            6.            6.            6.
(1,0,1,*)              7.            7.            7.            7.            7.
(1,1,0,*)              6.            6.            6.            6.            6.
(1,1,1,*)              7.            7.            7.            7.            7.
(1,2,0,*)              6.            6.            6.            6.            6.
(1,2,1,*)              7.            7.            7.            7.            7.

(0,0,0,*)              8.            9.           10.           11.           12.
(0,0,1,*)              8.            9.           10.           11.           12.
(0,1,0,*)              8.            9.           10.           11.           12.
(0,1,1,*)              8.            9.           10.           11.           12.
(0,2,0,*)              8.            9.           10.           11.           12.
(0,2,1,*)              8.            9.           10.           11.           12.
(1,0,0,*)              8.            9.           10.           11.           12.
(1,0,1,*)              8.            9.           10.           11.           12.
(1,1,0,*)              8.            9.           10.           11.           12.
(1,1,1,*)              8.            9.           10.           11.           12.
(1,2,0,*)              8.            9.           10.           11.           12.
(1,2,1,*)              8.            9.           10.           11.           12.

(0,0,0,*)             18.           19.           20.           21.           22.
(0,0,1,*)             19.           20.           21.           22.           23.
(0,1,0,*)             19.           20.           21.           22.           23.
(0,1,1,*)             20.           21.           22.           23.           24.
(0,2,0,*)             20.           21.           22.           23.           24.
(0,2,1,*)             21.           22.           23.           24.           25.
(1,0,0,*)             19.           20.           21.           22.           23.
(1,0,1,*)             20.           21.           22.           23.           24.
(1,1,0,*)             20.           21.           22.           23.           24.
(1,1,1,*)             21.           22.           23.           24.           25.
(1,2,0,*)             21.           22.           23.           24.           25.
(1,2,1,*)             22.           23.           24.           25.           26.

(0,0,*)            100.          105.
(0,1,*)            105.          110.
(0,2,*)            110.          115.
(1,0,*)            105.          110.
(1,1,*)            110.          115.
(1,2,*)            115.          120.

(0,*)            205.          215.          225.
(1,*)            215.          225.          235.

(0,*)            645.
(1,*)            675.

1320.

例子4 :在matlab中,純for循環速度最慢。而一半for循環+一半向量化的速度最快,Forcal中也是如此

!using["math","sys"];
mvar:
(:p1,p2,p3,a,b)=
{
    oo{
        a=array[1000].rand(),
        b=array[1000].rand(),
        p1=array[1000,1000],
        p2=array[1000,1000],
        p3=array[1000,1000],
        t0=clock(),
        ndgrid(a,b,&A,&B),
        p1.=A+B
    },
    printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},

    lena=FCDLen(a),
    lenb=FCDLen(b),
    t0=clock(),
    m = lenb-1, (m>=0).while{
        oo{p2(m,neg) = a+rn[b(m)]},
        m--
    },
    printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},

    t0=clock(),
    m = lenb-1, (m>=0).while{
        n = lena-1, (n>=0).while{
            //p3(m,n) = a(n)+b(m), //用這句還要慢一些
            A(p3,m,n) = A(a,n)+A(b,m),
            n--
        },
        m--
    },
    printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}
};

    結果

ndgrid: 3.2001e-002
for1: 1.4999e-002
for2: 1.86

例子5 :一段程序的Forcal實現:

//用C++代碼描述為:
s=0.0;
for(x=0.0;x<=1.0;x=x+0.0011)
{
   for(y=1.0;y<=2.0;y=y+0.0009)
   {
     s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   }
}

    (1)數組求和函數Sum

!using["math","sys"];
mvar:
t=clock(),
oo{
    ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn (2)))))),0]
};
[clock()-t]/1000;

    結果:

1008606.64947441
0.625   //時間

    (2)求和函數sum

f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
sum["f",0,1,0.0011,1,2,0.0009];

    結果:

1008606.64947441
0.719   
//時間

    (3)while循環

mvar:
t=sys::clock();
s=0,x=0,
while{x<=1,
//while循環算法;
    y=1,
    while{y<=2,
        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
        y=y+0.0009
    },
    x=x+0.0011
},
s;
[sys::clock()-t]/1000;

    結果:

1008606.64947441
0.734   //時間

例子6 :arrayfun應用:曲面擬合求參數

模型方程:z=107.6713*x^(b+0.3283)/(x^0.3283+24.3184*y^0.0093)

x y z數據如下:
5 0.5 24.791883
5 3.635 24.72700666
5 6.642 24.57905574
5 10.69 24.50235246
10 0.5 48.9584826
10 3.635 48.68938331
10 6.642 48.4001109
10 10.69 48.2125388
15 0.5 72.0477931
15 3.635 71.37663586
15 6.642 71.12836949
15 10.69 71.04465675
20 0.5 94.8989091
20 3.635 93.86373622
20 6.642 92.94844038
20 10.69 92.26085095

Forcal代碼:

!using["fcopt","math","io"];
f(x,y,z::b)=[107.6713*x^(b+0.3283)/(x^0.3283+24.3184*y^0.0093)-z]^2;
ff(bb ::x,y,z,b)= b=bb, oo{arrayfun[HFor("f"),x,y,z].Sum[]};
//目標函數定義
main(::x,y,z)=
{
    oo{ x=array[16].arrayns{"5, 5, 5, 5, 10, 10, 10, 10, 15, 15, 15, 15, 20, 20, 20, 20"},
        y=array[16].arrayns{"0.5 3.635 6.642 10.69 0.5 3.635 6.642 10.69 0.5 3.635 6.642 10.69 0.5 3.635 6.642 10.69"},
        z=array[16].arrayns{"24.791883 24.72700666 24.57905574 24.50235246 48.9584826 48.68938331 48.4001109 48.2125388 72.0477931 71.37663586 71.12836949 71.04465675 94.8989091 93.86373622 92.94844038 92.26085095"}
    },
    Opt[HFor("ff")]
//優化
};

結果:

0.7335293627980047 47.0803238072807

例子7 :arrayfun應用:向量元素乘

    有三個向量,分別是A=[1;2]; B=[3;4]; C=[5;6]; 需要對A,B,C中的元素進行乘法,需要考慮所有的可能性。這里1*3*5;1*3*6; 2*3*5; 2*3*6; 1*4*5; 1*4*6; 2*4*5; 2*4*6。

Forcal代碼:

!using["math","sys"];
mvar:
f(x,y,z)=x*y*z;
oo{
    A=array[2].SA[0:1,2],
    B=array[2].SA[0:3,4],
    C=array[2].SA[0:5,6],
    ndgrid[A,B,C,&x,&y,&z],
    arrayfun[HFor("f"),x,y,z].outm(),
//A,B,C中的元素乘法,考慮所有的可能性
    Sum[x*y*z,0]
//結果相加
};

結果:
            15.            18.            20.            24.            30.            36.            40.            48.

231.

例子8

 


版權所有© Forcal程序設計 2002-2010,保留所有權利
E-mail: [email protected]
  QQ:630715621
最近更新: 2011年02月05日

欧冠杯