2010年10月14日木曜日

// 水面（境界あり） (4)

コメント欄のプログラムをＳｃｉｌａｂで実行すると、下のような図が現れます。
（Ｓｃｉｌａｂ上では、なぜか横に傾いており、下の図は９０度傾けています。）
児島湖（のような形の池）を想定しています。
左の２つの流入河川から、池へ水が流れて行き、右上の出口から出て行きます。
水の流れは、本当は二次の微分方程式で表されるのですが、ここでは一次だけです
（単純に拡散しているだけ）。

※　今後、二次の式（ラプラシアンという）に、アップグレードしてゆきたいと思っています。

投稿者 時刻:

1 件のコメント:

  1. kf332010年10月14日 22:36

    // 水面（境界あり） (4)

    //流入と流出あり
    //三角池（児島湖）
    clear

    n=32
    K=32 //K 繰り返し回数


    Z(n,n)=0
    // Z(8,8)=100 //初期値
    Z1(n,n)=0
    Z2(n,n)=0

    fX(n-1,n-1)=0 //ベクトル表示の初期値
    fY(n-1,n-1)=0

    W=[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0
    0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0
    0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0
    0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0
    0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0
    0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0
    0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
    0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
    0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
    0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
    0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
    0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
    0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]




    // Z1(8-1:8+1,8-1:8+1)=W(8-1:8+1,8-1:8+1) .* Z(8,8)/sum(W(8-1:8+1,8-1:8+1))



    for K=1:K

    for i=3:n-2,
    for j=3:n-2,

    if sum(W(i-1:i+1,j-1:j+1))=0 then, Z1(:,:)=0, else,
    Z1(i-1:i+1,j-1:j+1)=W(i-1:i+1,j-1:j+1) .* Z(i,j)/sum(W(i-1:i+1,j-1:j+1)),end

    Z2=Z2+Z1,
    Z1(:,:)=0,

    end,
    end,

    Z=Z2,
    Z2(:,:)=0,

    Z(3,27:30)=0 //水の出口から流出
    Z(3:6,3)=10 //水の入り口がら流入
    Z(11:12,3)=10 //水の入り口がら流入

    //--------------------------プロット
    // XX=1:n,YY=1:n,
    // clf
    // plot3d1(XX,YY,Z),
    //-------------------



    XX=1:n-1
    YY=1:n-1

    for i=2 : n-2,
    for j=2 : n-2,
    fX(i,j)=( sum(Z(i-1:i+1, j-1 )) - sum(Z(i-1:i+1, j+1 )) )/3 ,
    fY(i,j)=( sum(Z(i-1 , j-1:j+1)) - sum(Z(i+1 , j-1:j+1)) )/3 ,
    end,
    end,

    clf,
    // xgrid(4)
    champ(XX,YY,fY,fX),

    fY(:,:)=0,
    fX(:,:)=0,


    end

    返信削除

登録: コメントの投稿 (Atom)