Return to rho CVS log

Up to [Research Unix] / researchv10dc / cmd / matlab / exec

Annotation of researchv10dc/cmd/matlab/exec/rho, revision 1.1

1.1     ! root        1: //Conductivity example.
        !             2: //Parameters ---
        !             3:    rho       //radius of cylindrical inclusion
        !             4:    n         //number of terms in solution
        !             5:    m         //number of boundary points
        !             6: //initialize operation counter
        !             7:    flop = <0 0>;
        !             8: //initialize variables
        !             9:    m1 = round(m/3);   //number of points on each straight edge
        !            10:    m2 = m - m1;       //number of points with Dirichlet conditions
        !            11:    pi = 4*atan(1);
        !            12: //generate points in Cartesian coordinates
        !            13:    //right hand edge
        !            14:    for i = 1:m1, x(i) = 1; y(i) = (1-rho)*(i-1)/(m1-1);
        !            15:    //top edge
        !            16:    for i = m2+1:m, x(i) = (1-rho)*(m-i)/(m-m2-1); y(i) = 1;
        !            17:    //circular edge
        !            18:    for i = m1+1:m2, t = pi/2*(i-m1)/(m2-m1+1); ...
        !            19:       x(i) = 1-rho*sin(t);  y(i) = 1-rho*cos(t);
        !            20: //convert to polar coordinates
        !            21:    for i = 1:m-1, th(i) = atan(y(i)/x(i));  ...
        !            22:       r(i) = sqrt(x(i)**2+y(i)**2);
        !            23:    th(m) = pi/2;  r(m) = 1;
        !            24: //generate matrix
        !            25:    //Dirichlet conditions
        !            26:    for i = 1:m2, for j = 1:n, k = 2*j-1; ...
        !            27:       a(i,j) = r(i)**k*cos(k*th(i));
        !            28:    //Neumann conditions
        !            29:    for i = m2+1:m, for j = 1:n, k = 2*j-1; ...
        !            30:       a(i,j) = k*r(i)**(k-1)*sin((k-1)*th(i));
        !            31: //generate right hand side
        !            32:    for i = 1:m2, b(i) = 1;
        !            33:    for i = m2+1:m, b(i) = 0;
        !            34: //solve for coefficients
        !            35:    c = A$b
        !            36: //compute effective conductivity
        !            37:    c(2:2:n) = -c(2:2:n)
        !            38:    sigma = sum(c)
        !            39: //output total operation count
        !            40:    ops = flop(2)

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.