วันศุกร์ที่ 15 พฤษภาคม พ.ศ. 2563

บทที่ 5 การทำงานซ้ำกับการประยุกต์ในปัญหาเมมเบรน


บทที่ 5

วิธีทำซ้ำประยุกต์ใช้กับปัญหาเมมเบรน

5.1 ส่วนนำ
            ในเอกสารเสริมนี้จะได้สืบเสาะหาคำตอบจากปัญหาเมทริกซ์ซึ่งพบโดยทั่วไปในการจัดทำโมเดลของเมมเบรน ความสมดุล การแจกแจงของอุณหภูมิ การไหลของของไหลผ่านเยื่อรูพรุน( porous media) .
ในครั้งแรกจะกำหนดปัญหาและแก้ปัญหาโดยใช้การแก้สมการเชิงเส้นแบบสปาร์สของไซแลบ  จากนั้นจะสืบเสาะความหมายการเพิ่มความเร็ว โดยเฉพาะอย่างยิ่งจะเห็นว่าวิธีการทำงานซ้ำ  วิธีการคอนจูเกต
กราเดียน (conjugate gradient method) สำหรับปัญหาการแก้ปัญหาขนาดใหญ่   เอกสารเสริมนี้จะให้ประสบการณ์แก่ผู้ใช้ในการพล็อตกราฟของฟังก์ชัน 2 ตัวแปร (surface graphs).
            มีรหัสคำสั่งจำนวนมากในเอกสารเสริมนี้ ซึ่งนำวิธีการที่ละเอียดมาใช้  จึงขอแนะนำ การ cut และการ paste นำรหัสคำสั่งเข้ามาใช้ในไซแลบและรันให้ทำงานตามคำสั่งที่ระบุในรหัสคำสั่ง

5.2 ติดตั้งเมทริกซ์
            ครั้งแรกสุดตั้งปัญหาหนึ่งซึ่งรูปแบบของเมมเบรนหนึ่งด้วยการการโหลด (loading) เช่นโหลดมวล (mass loads) หรือ ความดัน  เมมเบรนมีรูปแบบกริด 2 มิติโดยประมาณ
ด้วยการติดตั้งกริดดังนี้

n = 10; m = 15;
x = 0:1/(n-1):1;
y = 0:1/(m-1):1;
u = x’*y;

            ผู้ใช้สามารถพล็อตค่า u โดยใช้คำสั่งการพล็อตดังนี้

plot3d(x,y,u)

หมุนภาพการพล็อต (โดยใช้ป่มการหมุนของภาพในวินโดว์แสดงกราฟ) เพื่อจะให้ได้แนวคิดที่ดีของรูปร่างของพื้นผิว  ผู้ใช้อาจต้องการทดลองใช้การพล็อตรูปแบบที่แตกต่างเช่น  shaded plot โดยใช้คำสั่ง  plot3d1 ดังนี้
plot3d1(x,y,u)

5.3 การติดตั้งค่าขอบ
            ความต้องการที่สร้างโมเดลในสถานะการณ์ซึ่งเมมเบรนกำหนดตายตัวรอบขอบ  ต่อไปจะได้ระบุความสูงของโหนดขอบ (edge nodes) โดยใช้ฟังก์ชันหนึ่ง  ในขณะแรกกำหนดขอบเขตที่จะให้ได้รูปร่างอานม้ากำหนดโดยฟังก์ชัน ดังนี้

function u = saddle_function(x,y)
     //
    // Saddle point function to be used as
    // boundary and exact solution
    //
    u = (x-0.5)^2 - (y-0.5)^2
endfunction

           
ให้โหลดฟังก์ชันนี้เข้าไปในไซแลบ โหนดรอยต่อหรือเขตแดนที่ตรงกับ ค่า  i, j  เมื่อ i = 1, i = n,   j = 1 และ j = m. ดังนั้นให้กำหนดโหนดเหล่านั้นเพื่อให้มีค่ากำหนดโดยฟังก์ชัน

saddle function.
for i=1:n
   for j=1:m
      // If statement to determine boundary nodes
       If( i==1 | i==n | j==1 | j==m) then
       u(i,j) = saddle_function(x(i),y(j));
    else
     u(i,j) = 1.0;
   end
  end
end

จากนี้ให้ทดลองพล็อตกริดฟังก์ชัน u  การพล็อตต่อไปใช้คำสั่ง  plot3d


5.4 ประยุกต์กับแรงตึง (Tension)
            เมื่อต้องการประยุกใช้แรงตึงกับเมมเบรน ทำได้โดยการกำหนดดังนี้
fij = nm/4 [ui+1j + ui1j + uij+1 + uij1] uij
            จากนี้ คือ fij แทนค่าแรงทางแนวตั้ง ประยุกต์ใช้กับโหนด  i, j. แรงทางแนวตั้งนั้นโดยทั่วไปมักตีความเป็นแรงความโน้มถ่วงหรือแรงจากความดัน ดังนั้นทำให้มีระบบของสมการเชิงเส้น
ui+1j + ui1j 4uij + uij+1 + uij1 = 4fij /nm
            จำนวนของโหนดคือ nm,และดังนั้นเมทริกซ์ที่สอดคล้องกันก์คือ nm โดย nm.
.ให้กำหนดเมทริกซ์สปาร์สซึ่งครอบคลุมสมการเหล่านี้ทั้งหมด สำหรับโหนดขอบจะติดตั้งแถวที่สอดคล้องตรงกันของเมทริกซ์เพื่อให้มีเฉพาะตามแนวทะแยง เท่ากับหนึ่ง และกำหนดทางขวาเท่ากับค่าขอบหรือเขตแดน  โดยเริ่มต้นเราจำเป็นต้องกำหนดสปาร์สเมทริกซ์ โดยใช้ฟังก์ชันการสร้างเมทริกซ์ต่อไปนี้
//--------------------------------------------
// Define Functions
//--------------------------------------------
function [A,b] = create_matrix(n,m)
//          
// Create a matrix associated with a grid
// nodes with x and y coordinates given as
// input, connected by springs with nodes
// on the boundary set the function boundary
//
x = 0:1/(n-1):1;
y = 0:1/(m-1):1;
A = speye(n*m,n*m);
b = zeros(n*m,1);
for i=1:n
   for j=1:m
     index = i + (j-1)*n;
     if( i==1 | i==n | j==1 | j==m) then // Boundary nodes
        b(index) = saddle_function(x(i),y(j));
        A(index,index) = 1.0;
      else // Interior Nodes
         b(index) = 4*force/(n*m);
         A(index,index) = -4.0;
         A(index,index+1) = 1.0;
         A(index,index-1) = 1.0;
         A(index,index+n) = 1.0;
         A(index,index-n) = 1.0;
       end
     end
   end
endfunction

            ให้โหลดฟังก์ชันนี้เข้าไปในไซแลบ และสร้างเมทริกซ์ A และทางขวามือ b ตรงกับ a5 โดย 5 grid โดยการใช้คำสั่ง
            [A,b]=create_matrix(5,5)
กรณีนี้จะก่อให้เกิดข้อผิดพลาดเนื่องจากฟังก์ชันเมทริกซ์ที่สร้างขึ้นต้องการชุดตัวแปรบังคับ  โดยกำหนดให้เป็น 0 และทดลองใช้คำสั่ง
force = 0.0;
[A,b]=create_matrix(5,5)

            เมทริกซ์ A มีใหญ่อย่างไร และมีโครงสร้างใด  ต่อไปเป็นรหัสคำสั่งสำหรับฟังก์ชัน spy ซึ่งมีประโยชน์เพื่อมองภาพโครงสร้างของเมทริกซ์

function spy(A)
// Mimic the Matlab spy command
// Draw the non zero components of a matrix
[i,j] = find(A~=0)
[N,M] = size(A)
// wrect=[x,y,w,h] frect=[xmin,ymin,xmax,ymax]
xsetech([0,0,1,1],[1,0,M+1,N])
xrects([j;N-i+1;ones(i);ones(i)],ones(i));
// [x,y,w,h]
xrect(1,N,M,N);
endfunction

5.5 การแก้ปัญหา
            จากนี้มีแนวทางในการสร้างเมทริกซ์ และด้านขวามือ  ดังนั้นตอนนี้ก็เพียงแต่แก้สมการ ดังแสดงรหัสคำสั่งที่คิดขึ้นเพื่อแก้สมาการ โดยได้เพิ่มรหัสคำสั่งเพิ่มเติม  ซึ่งในครั้งใดที่รหัสคำสั้งได้เพิ่มส่วนรัสพิเศษ ซึ่งเวลาที่รหัสคำสั่งได้กำหนด  wrapper สำหรับผู้แก้ปัญหา ดังนั้นเราสามารถเปลี่ยนแปลงตัวแก้ปัญหาต่อมาในเอกสารเสริมนี้.

function [U] = solve_equation(A,b)
//
// Wrapper for linear solver
//
U = A\b
endfunction
function run_problem(n,m)
//
Computational Science, Scilab Tutorials 41
// Setup a problem, create matrix,
// solve and then plot result.
//
// Input:
// n,m, number of grid points in x and y direction.
//
x = 0:1/(n-1):1;
y = 0:1/(m-1):1;
printf(’Create Matrix:’)
timer()
[A,b] = create_matrix(n,m);
printf(’ time = %g\n’,timer())
printf(’Solve Matrix problem:’)
timer()
U = solve_equation(A,b)
printf(’ time = %g\n’,timer())
u = matrix(U,n,m);
// Plot the solution
xset("colormap",jetcolormap(64))
xbasc(); xselect();
plot3d1(x,y,u)
endfunction

            จากนี้ควรจะหมดรหัสคำสั่งแล้ว ซึ่งก็เพียงแต่ต้องรันฟังก์ชันที่รันปัญหา  ตัวอย่างเช่นเพื่อที่จะ รันรหัสคำสั่งสำหรับ10 โดย 15 กริด ใช้คำสั่ง  run_problem(10,15) เมื่อรันส่วนของรหัสคำสั่งที่แล้ว ทำให้ได้มาซึ่ง การพล็อตต่อไปนี้
            นี่ทำได้ดีเท่าที่สามารถเห็นได้ว่า เมมเบรนกำหนดตนายตัวไปยังขอบ และภายในเมมเบรนดูเหมือนว่าดึงออกไปอย่างเท่าเทียมกันในแต่ละทิศทาง
การบันทึกการสร้างการพล็อตเมมเบรนด้วยเงื่อนไขขอบ saddle เดียวกัน แต่ด้วยประยุกตแรง 5  
การพล็อตกราฟเมมเบรนด้วยเงื่อนไขขอบเป็นศูนย์  และประยุกต์ให้แรง 5  นั้นจำเป็นต้องเปลี่ยนการสร้างฟังก์ชันเมทริกซื.
            จากนี้ทดลองการใช้กริดที่โตขึ้น  ทดลองแก้ปัญหากล่าวคตือกริด 40 by 40 หรือ 50 by 50.
การบันทึกต่อไปว่าใช้เวลานานเท่าใดที่จะสร้างเมทริกซ์ที่เกี่ยวพันธ์กับกริด  50 by50 grid.  ใช้เวลานานเท่าใดในการแก้สมการที่เกี่ยวพันธ์นั้น  และเมทริกซ์มีขนาดโตเท่าใด.

5.6 การสร้างเมทริกซ์ที่เร็วกว่า
            ที่ผ่านมาใช้เวลาในการสร้างเมทริกซ์นานกว่าการแก้สมการเมทริกซ์ คิดได้ว่าอาจกำลังทำอะไรบางอย่างผิดพลาด  ในแต่ละส่วนของเมทริกซ์ (sparse matrix) ได้ปรับปรุงไปคือการใช้เวลามาก ในการนำไปใช้ต่อไปของฟังก์ชันการสร้างเมทริกซ์  โดยใช้ฟังก์ชัน diag(sparse) เพื่อติดตั้งเมทริกซ์ ให้โหลดการนำไปใช้ใหม่ต่อไปนี้ในฟังก์ชันการสร้างเมทริกซ์เข้าไปในไซแลบ

function [A,b] = create_matrix(n,m)
//
// Create a matrix associated with a grid
// nodes with x and y coordinates given as
// input, connected by springs with nodes
Computational Science, Scilab Tutorials 43
// on the boundary set the function boundary
//
x = 0:1/(n-1):1
y = 0:1/(m-1):1
//-----------------------
// Setup Matrix
//----------------------
A = spzeros(n*m,n*m);
//
// Setup interior problem
//
d0 = -4*ones(n*m,1);
dmn = ones((m-1)*n,1)
dpn = ones((m-1)*n,1)
dm1 = ones(n*m-1,1)
dp1 = ones(n*m-1,1)
//
// Apply Boundary condition
//
d0(n:n:n*m) = 1
d0(1:n:n*m) = 1
d0(1:n) = 1;
d0(n*(m-1):n*m) = 1
dmn(1:n:n*(m-1)) = 0
dmn(n:n:n*(m-1)) = 0
dmn(n*(m-2):n*(m-1)) = 0
dpn(1:n:n*(m-1)) = 0
dpn(n:n:n*(m-1)) = 0
dpn(1:n) = 0
dm1(n-1:n:n*m-1) = 0
dm1(n:n:n*m-1) = 0
dm1(1:n-1) = 0
dm1(n*(m-1):n*m-1) = 0
dp1(n+1:n:n*m-1) = 0
dp1(n:n:n*m-1) = 0
Computational Science, Scilab Tutorials 44
dp1(1:n-1) = 0
dp1(n*(m-1):n*m-1) = 0
A = diag(sparse(d0)) + diag(sparse(dm1),-1) + diag(sparse(dp1),1) + ...
diag(sparse(dmn),-n) + diag(sparse(dpn),n)
//----------------
// Setup rhs
//----------------
b = 4*force/(n*m)*ones(n*m,1);
//
// Apply boundary condition
//
b(1:n) = feval(x,y(1),saddle_function);
b(n*(m-1)+1:n*m) = feval(x,y(m),saddle_function);
b(1:n:n*(m-1)+1) = feval(x(1),y,saddle_function)’;
b(n:n:n*m) = feval(x(n),y,saddle_function)’;
endfunction

ให้เชื่อมั่นในตัวเองว่ารหัสคำสั่งนี้ทำงานเหมือนกับการนำไปใช้ในคราวที่แล้ว.ให้ดำเนินการปัญหาของกริด 50 by 50   ให้เปรียบเทียบอัตราเร็วของการนำไปใช้ใหม่กับการนำไปใช้ที่แล้วมา รหัสคำสั่งที่ใช้ทำงานได้เร็วขึ้นหรือไม่ สามารถเห็นได้ว่าผลเฉลยตอนนี้ของระบบเชิงเส้นนั้นดีกว่าเวลาในการทำเนินการ

5.7 วิธีการแกรเดียนคอนจูเกต
            ส่วนนำเข้าสู่วิธีการทำงานซ้ำๆ ในส่วนนี้วิธีการคอนจูเกทซ์ทำงานได้ดีกับเมทริกซ์เช่นเมทริกซ์ A ซึ่งจริงแล้วไม่ได้มีความสมมาตร แต่ตราบเท่าที่แน่ใจว่าการคาดคะเนสอดคล้องกับเงื่อนไขขอบ แล้วโดยหลักการแล้วเมทริกซ์มีความสมมาตรเมื่อเทียบกับจุดกริดภายใน  วิธีกราเดียนคอนจูเกทไม่มีให้ใช้ในไซแลบโดยปริยาย แต่เป็นอัลกอริธึมที่ง่ายและสามารถนำมาใช้ผ่านทางรหัสคำสั่งดังต่อไปนี้ (นำมาจาก An Introduction to the Conjugate Gradient Method Without the Agonizing Pain”, by Jonathan Richard Shewchuk)
function [x,i] = conjugate_gradient(A,b,x0,imax,tol,iprint)
//
// Try to solve linear equation Ax = b using
// conjugate gradient method
//
// Input
// A: matrix or function which applies a matrix, assumed symmetric
// b: right hand side
// x0: inital guess
// imax: max number of iterations
// tol: tolerance used for residual
//
// Output
// x: approximate solution
// i: number of iterations
//
//------------------------------
// Check and set input arguments
//------------------------------
nargin = argn(2)
if nargin<2
error(’need at least 2 input arguments’)
end
if nargin<3
x0 = zeros(b);
end
if nargin<4
imax = 2000;
end
if nargin<5
tol = 1e-8
end
if nargin<6
iprint = 0
end
function y = apply_A(x)
if (typeof(A)==’function’) then
y = A(x)
else
y = A*x
end
endfunction
//-------------------------------
i=0
x = x0
r = b - apply_A(x)
d = r
rTr = r’*r
rTr0 = rTr
while (i<imax & rTr>tol^2*rTr0),
q = apply_A(d)
alpha = rTr/(d’*q)
x = x + alpha*d
if modulo(i,50)==0 then
r = b - apply_A(x)
else
r = r - alpha*q
end
rTrOld = rTr
rTr = r’*r
bt = rTr/rTrOld
d = r + bt*d
i = i+1
if modulo(i,iprint)==0 then
printf(’i = %g rTr = %20.15e \n’,i,rTr )
end
end
if i==imax then
warning(’max number of iterations attained’)
end
endfunction
Load this function into Scilab and try on a small matrix and rhs such as
n=21;
B = diag(sparse(ones(n-1,1)),-1) - 2*diag(sparse(ones(n,1))) + ...
diag(sparse(ones(n-1,1)),1);
c = ones(n,1);
by using the command
x = conjugate_gradient(B,c)

            เป็นการดีแต่สามารถเห็นรายละเอียดได้มากขึ้นถ้าใช้อาร์กิวเมนต์มากขึ้น
x = conjugate_gradient(B,c,zeros(c),50,1e-9,1)
            อาร์กิวเมนต์ที่ 3 จัดเป็นการคาดคะเนเริ่มต้นสำหรับแผนการทำงานซ้ำ (iterative scheme) อาร์กิวเมนต์ที่ 4 จัดเป็นขีดจำกัดของของจำนวนการทำงานซ้ำ อาร์กิวเมนต์ที่ 5 การคงทนต่อการจุดความสัมพันธ์ของส่วนเหลือ และอาร์กิวเมนต์ที่ 6 ให้เอ้าพุตการวินิจฉัยตามปกติ ผลลัพธ์แสดงให้เห็นว่ากำหลังสองของส่วนเหลือ r*r ลดลงอย่างช้าๆ จนกระทั้งการทำงานซ้ำที่ 11 ที่ลดลงไปเป็น0  อันเป็นวิธีการปกติ  นั้นสุดท้ายสามารถที่จะใช้วิธีแกรเดี๊อนคอนจูเกตเพื่อพยายามในการแก้ปัญหาเมมเบรน  การรันฟังก์ชันปัญหาเรียกฟังก์ชันที่แก้สมการ การใช้วีแกรเดียนคอนจูเกตจะเป็นต้องนำสมการที่แก้ปัญหามาใช้ใหม่.

function [U] = solve_equation(A,b)
//
// Get good initial guess. well at least one that
// satisfies BC
//
x0 = b
U = conjugate_gradient(A,b,x0,length(x0),1e-9,0)
Endfunction

การบันทึกการรันปัญหาเมมเบรนสำหรับกริด 50 by 50 grid, โดยใช้ฟังก์ชันการแก้สมการ ให้ข้อคิดเห็นต่ออัตราเร็วสัมพันธ์ของการนำมาใช้ใหม่ของรหัสคำสั่ง  ว่าดีกว่าการนำไปใช้ที่ผ่านมาหรือไม่ หาโดยประมาณว่าปัญหาใหญ่โตเท่าใดที่สามารถแก้ได้ภายใน 20 วินาทีด้วยคอมพิวเตอร์ด้วยการนำมาใช้ท้ายสุด.


5.8 การคูณเว็คเตอร์เมทริกซ์ที่เร็วกว่า
            ตามที่เป็นไปตามปกตินั้นไม่ใช้จะจบเรื่องกันไป สิ่งแรกสามารถที่จะนำวิธีการที่ดีกว่าไปใช้สำหรับการคำนวณเว็คเตอร์เมทริกซ์(ส่วนที่ยากที่สุดของอัลกอริธึมแกเดียนคอนจูเกต).  ในความเป็นจริงถ้าใช้วิธีการทำซ้ำเช่นวิธีแกรเดียนคอนจูเกต ก็เพียงแต่จำต้องจัดหาขันตอนวิธีสำหรับการคำนวณการคูณเว็คเตอร์เมทริกซ์  ดังนั้นต่อไปนี้ยังมีการนำขั้นตอนวิธีอื่นในการสร้างเมทริกซ์ไปใช้

function [A,b] = create_matrix(n,m)
//
// Create a matrix operator associated with a grid
// nodes with x and y coordinates given as
// input, connected by springs with nodes
// on the boundary set the function boundary
//
x = 0:1/(n-1):1
y = 0:1/(m-1):1
//-----------------------
// Setup Matrix
// A function is being created
// which will be returned
// instead of a matrix
//----------------------
function y = A(x)
//
// Setup finite difference operator for Laplacian
//
x = matrix(x,n,m)
y = x
J = 2:m-1
I = 2:n-1
y(I,J) = -4.0*x(I,J) + x(I+1,J) + x(I-1,J) + x(I,J+1) + x(I,J-1)
y = y(:)
endfunction
//----------------
// Setup rhs
//----------------
b = 4*force/(n*m)*ones(n*m,1);
//
// Apply boundary condition
//
b(1:n) = feval(x,y(1),saddle_function);
b(n*(m-1)+1:n*m) = feval(x,y(m),saddle_function);
b(1:n:n*(m-1)+1) = feval(x(1),y,saddle_function)’;
b(n:n:n*m) = feval(x(n),y,saddle_function)’;
Computational Science, Scilab Tutorials 49
Endfunction

            ในกรณีนี้ขั้นตอนวิธีคืนกลับไม่ใช่สปาร์สเมทริกซ์ A แต่เป็นขั้นตอนวิธีหนึ้งสำหรับประยุกต์ใช้ การบันทึกให้มองไปที่รหัสคำสั่งใหม่สำหรับขั้นตอนวิธีสร้างเมทริกซ์  อธิบายว่าขั้นตอนวิธีใดกำลังทำงาน โดยเฉพาะอย่างยิ่ง อธิบายเส้นต่อไปนี้

J = 2:m-1
I = 2:n-1
y(I,J) = -4.0*x(I,J) + x(I+1,J) + x(I-1,J) + x(I,J+1) + x(I,J-1)

การบันทึกการทดลองวิธีการใหม่นี้  ว่าปัญหาใหญ่เท่าใดที่สามารถแก้ได้ภายใน 20 วินาที

5.9 สรุป
            แม้ว่าตอนนี้มีที่ว่างสำหรับการปรับปรุง วิธีการทำงานซ้ำสามารถที่จะปรับปรุงให้ดีขึ้นได้ โดยใช้สิ่งที่เรียกว่าเงื่อนไขก่อนล่วงหน้า ด้วยการปรับปรุงเหล่านี้ เป็นไปได้ที่จะแก้ปัญหากริดขนาด  1000 by 1000 (1,000,000 ตัวไม่ทราบค่า)  เป็นเวลา 10 วินาทีบนเครื่องพีซีธรรมดาร
            อย่างที่กล่าวมาอาจจะตั้งคำถามว่าทำไมจำเป็นต้องแก้ปัญหากริดขนาด 1000 by 1000  ตัวอย่างเช่นในสถานะการณ์ทางธรณีกายภาพที่มักมีกรณีว่าที่ใดที่ต้องแก้ปํญหาการไหลไม่กี่เมตร (โดยปริยาย) บนพื้นหลายกิโลเมตร นั่นคือช่วงของมาตรมากกว่าพัน

ไม่มีความคิดเห็น:

แสดงความคิดเห็น