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

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.