โปรแกรมช่วยหาคำตอบระบบสมการไม่เชิงเส้นด้วย Newton-Raphson Method

สวัสดีครับ ในบทความนี้ จะขอนำเสนอการประยุกต์คอนสตรัคเตอร์แบบที่ใช้การโหลดสมาชิกของเมตริกซ์ จากไดนามิกส์อาเรย์  โดยจะนำเสนอผ่านระเบียบวิธีการหาคำตอบของระบบสมการไม่เชิงเส้นด้วย Newton-Raphson Method ผมจะขออธิบายแบบสั้นๆก่อนนำเข้าสู่โค้ดโปรแกรมภาษา C++ โดยปกติระบบสมการเชิงเส้น เราสามารถเขียนอยู่ในรูปเมตริกซ์ได้เป็น

[A][x] = [b]

โดย เมตริกซ์ A เป็นเมตริกซ์ที่มีสมาชิกของเมตริกซ์เป็นค่าคงที่ ดังแสดงในรูปที่ 1

ซึ่งเมตริกซ์ A จะเป็นเมตริกซ์ที่คงที่ ดังนั้นจึงสามารถใช้คลาส CMatrix อ่านข้อมูลสมาชิกของเมตริกซ์ A และ b ดังแสดงในบทความ การหาคำตอบของระบบสมการเชิงเส้นได้ แต่ในกรณีเป็นระบบสมการไม่เชิงเส้น (Nonlinear Equation System) สมาชิกของเมตริกซ์ A จะเป็นฟังก์ชันของตัวแปรที่ต้องการหาคำตอบ ดังแสดงในรูปที่ 2 การหาคำตอบของระบบสมการไม่เชิงเส้นจึงต้องใช้การหาคำตอบของระบบสมการด้วย Newton-Raphson Method 

Newton-Raphson Method จะมีขั้นตอนดำเนินการดังนี้

1. เขียนสมการทั้งหมดให้อยู่ในรูปแบบดังนี้

จากตัวอย่าง 3 สมการจะสามารถเขียน f1 , f2 , f3 ได้ดังนี้

2. การหาคำตอบของระบบสมการ จะเริ่มจากการเขียนระบบสมการในรูปแบบของอนุกรมเทย์เลอร์ โดยสามารถเขียนในรูปเมตริกซ์ได้ดังนี้

               

เมตริกซ์ J จะเรียกว่า Jacobian matrix ซึ่งสมาชิกต่างๆหาได้ดังนี้

จากตัวอย่าง 3 สมการ จะเขียนระบบสมการได้ดังนี้

3. ทำการหาคำตอบของระบบสมการนี้ โดยการกำหนดค่าเริ่มต้น  

แทนลงในเมตริกซ์ J และ f และหาคำตอบของระบบสมการในข้อ 2 ดังกล่าว

4. คำนวณหาผลลัพธ์ใหม่จาก  

5. ตรวจสอบผลลัพธ์โดยการแทนคำตอบจากข้อ 4 ลงใน เมตริกซ์ f และตรวจสอบค่าของสมาชิกในเมตริกซ์ โดยหากคำตอบที่แทนลงในเมตริกซ์ f ทำให้สมาชิกในเมตริกซ์ f มีค่าเท่ากับ 0  จะได้ว่า คำตอบดังกล่าวเป็นคำตอบของระบบสมการ แต่ในทางปฏิบัติจะกำหนดเงื่อนไขการลู่เข้าหาคำตอบไว้ โดยหากค่าในเมตริกซ์ f ทุกค่าต่ำกว่าที่กำหนดจะถือว่าชุดคำตอบนี้เป็นคำตอบของระบบสมการ แต่หากไม่ใช่จะย้อนกลับไปหาคำตอบของระบบสมการในข้อ 2 ต่อไปและจะทำซ้ำจนกระทั่งลู่เข้าสู่คำตอบของระบบสมการ

จากที่ได้อธิบายมาคร่าวๆ เราจะเห็นจุดที่จะนำเอาสมาชิกเมตริกซ์จากอาเรย์ 2D เข้ามาประยุกใช้ได้ นั่นคือขั้นตอนที่ 2 จากนั้นกำหนดคำตอบเริ่มต้นในข้อ 3 และทำการหาคำตอบของระบบสมการตามข้อ 2 โดยใช้ฟังก์ชัน SolveLU หรือ SolveGauss ได้ เรามาดู โค้ดโปรแกรมกันเลยครับ

#include "stdafx.h"

#include "Matrix.h"

#include "math.h"

int _tmain(int argc, _TCHAR* argv[])

{

        int n ;

int count ;

int row , col ;

double **J ;

double **f ;

double **x ;

double tol = 0.000001 ;

row = 3 ;

col = 3 ;

count = 0 ;

J = new double*[row];

f = new double*[row];

x = new double*[row];

for(int i = 0 ; i < row ; i++)

{

J[i] = new double[col] ;

f[i] = new double[col] ;

x[i] = new double[col] ;

}

//เดาคำตอบเริ่มต้น

x[0][0] = 1 ;   //x

x[1][0] = 1 ; //y

x[2][0] = 1 ; //z

       //คำนวณสมาชิกเมตริกซ f เริ่มต้น

f[0][0] = -(4*x[0][0]+x[1][0]*x[1][0] +x[2][0] -11) ;  //-f1:-(4x+y^2+z-11)

f[1][0] = -(x[0][0]+4*x[1][0] +x[2][0]*x[2][0] -18) ;  //-f2:-(x+4y+z^2-18)

f[2][0] = -(x[0][0]*x[0][0]+x[1][0] +4*x[2][0] -15) ;  //-f3:-(x^2+y+4z-15)

do

{      

               count++;

       //คำนวณเมตริซ์ J

      J[0][0] = 4 ;

      J[0][1] = 2*x[1][0] ;  //2y

      J[0][2] = 1 ;

      J[1][0] = 1 ;

      J[1][1] = 4 ;

      J[1][2] = 2*x[2][0] ; //2z

      J[2][0] = 2*x[0][0] ; //2x

      J[2][1] = 1 ;

      J[2][2] = 4 ;

     CMatrix A(3,3,J) ;

     CMatrix b(3,1,f) ;

     CMatrix dx ;                      //สร้าง object dx ไว้สำหรับเก็บคำตอบที่ Solve ได้

     dx.SolveLU(&A,&b);       // Solve ระบบสมการและเก็บคำตอบไว้ที่เมตริกซ dx

       

             // กำหนดคำตอบที่ได้ลงในเมตริกซ์ x 

     x[0][0] += dx.GetData(0,0);

     x[1][0] += dx.GetData(1,0);

     x[2][0] += dx.GetData(2,0);

             //คำนวณสมาชิกเมตริกซ f ใหม่

     f[0][0] = -(4*x[0][0]+x[1][0]*x[1][0] +x[2][0] -11) ;  //-f1:-(4x+y^2+z-11)

     f[1][0] = -(x[0][0]+4*x[1][0] +x[2][0]*x[2][0] -18) ;  //-f2:-(x+4y+z^2-18)

     f[2][0] = -(x[0][0]*x[0][0]+x[1][0] +4*x[2][0] -15) ;  //-f3:-(x^2+y+4z-15)     

}

        // check เงื่อนไขการลู่เข้าของคำตอบ

while (fabs(f[0][0]) > tol || fabs(f[1][0]) > tol || fabs(f[2][0]) > tol);   

printf("inter %d\tsolution are :\n",count);

printf("x = %.3f\n",x[0][0]);

printf("y = %.3f\n",x[1][0]);

printf("z = %.3f\n\n",x[2][0]);

printf("f1 = %.7f\n",f[0][0]);

printf("f2 = %.7f\n",f[1][0]);

printf("f3 = %.7f\n\n",f[2][0]);

// คืน memory ให้ระบบ

for (int i = 0 ; i < row ; i++)

{

delete [] J[i] ;

delete [] f[i] ;

delete [] x[i] ;

}

delete [] J ;

delete [] f ;

delete [] x ;

J = NULL ;

f = NULL ;

x = NULL ;

        

return 0 ;

}

ผลการประมวลผลโปรแกรมภาษา c++ ดังกล่าวสำหรับใช้หาคำตอบระบบสมการที่ไม่เป็นเชิงเส้นตามตัวอย่างแสดงได้ดังรูป

คำตอบของระบบสมการ

จากตัวอย่างโค้ดจะเห็นว่าในส่วนของการสร้างเมตริกซ์ J และ f เป็นการประยุกต์ใช้วิธีการประกาศคอนสตรัคเตอร์ตามที่ได้นำเสนอในบทความที่ผ่านมาครับ ซึ่งก็จะช่วยให้วิศวกรสามารถแก้ปัญหาในระบบต่างๆได้โดยง่าย

         อยากจะกล่าวให้เห็นปัญหาของการหาคำตอบของระบบสมการที่มากกว่า 3 สมการ ซึ่งท่านต้องสร้างเมตริกซ์ J และ f มากขึ้น ซึ่งการสร้างเมตริกซ์ตามที่ผมได้นำเสนออาจจะทำได้ช้าและยุ่งยาก ซึ่งเราอาจจะประยุกต์การสร้างฟังก์ชันสำหรับคำนวณ สมาชิกในเมตริกซ์ขึ้นมาซึ่งทำได้รวดเร็วและเป็นแบบแผน ทำให้เราสามารถประยุกต์ใช้โค้ดโปรแกรมเดียวกันนี้กับระบบสมการหลายๆระบบได้โดยการเปลี่ยนแปลงโค้ดเพียงเล็กน้อย ซึ่งจะข้อนำเสนอในบทความถัดๆไปนะครับ เพื่อให้เห็นประโยชน์ของการเขียนโปรแกรมแบบ OOP 

ตัวอย่างการโหลดสมาชิกเมตริกซ์จากอาเรย์ 2D มาใช้ใน Class CMatrix ในโปรแกรมภาษา C++

สวัสดีครับ บทความนี้ผมจะขอนำเสนอการโหลดสมาชิกของเมตริกซ์จากอาเรย์ 2D มาใช้ในคลาส CMatrix กันครับ เริ่มจาก จากบทความแนะนำ คลาส CMatrix  ได้อธิบาย คอนสตรัคเตอร์ตัวหนึ่งที่ใช้ในการโหลดสมาชิกของเมตริกซ์จากอาเรย์ 2D ซึ่งมีรูปแบบการประกาศ Object ดังนี้ครับ

CMatrix A(const int r , const int c , double **pD)

มาพิจารณาพารามิเตอร์ที่ต้องส่งผ่านคอนสตรัคเตอร์ตัวนี้ ซึ่งจะประกอบไปด้วย

1. จำนวนแถวของเมตริซ์ r

2. จำนวนหลักของเมตริกซ์ c

3. อาเรย์ 2D ซึ่งต้องเป็นแบบ pointer 

สมมุติตัวอย่าง ต้องการโหลดสมาชิกเมตริกซ์ขนาด 3x3 จากตัวแปรอาเรย์ชื่อ A เราจะเขียน code โปรแกรมภาษา C++ ในส่วนนี้ได้ดังนี้

        int row , col ;

double** data ;

row = 3 ;

col = 3 ;

data = new double*[row];

for (int i = 0 ; i < row ; i++)

{

data[i] = new double[col] ;

}

for (int j = 0 ; j < row ; j++)

{

for(int k = 0 ; k < col ; k++)

data[j][k] = j+k ;

}

cout<<"Matrix A"<<endl;

CMatrix A(row,col,data);

A.Show();

cout<<endl ;

cout<<"Matrix B"<<endl;

CMatrix B(3,3,data) ;

B.Show();

ผลลัพธ์การทำงานใน code ส่วนนี้จะแสดงได้ดังนี้

ผลการทำงานของโค้ดโปรแกรม

จากตัวอย่างโค้ดโปรแกรมภาษา C++ ที่ประมวลผลเมตริกซ์นี้ ผมได้นำเสนอ 2 รูปแบบ คือ ท่านสามารถผ่านขนาดของเมตริกซ์ด้วยตัวแปรก็ได้ หรือจะผ่านด้วยค่าก็ได้ เนื่องจาก คอนสตรัคเตอร์นี้ กำหนดชนิดของตัวแปรขนาดเป็นแบบ คงที่ (const int r , const int c , double **pD) เอาหล่ะครับ จบตรงนี้เราก็ได้ object ของคลาส CMatrix ที่พร้อมจะนำไปใช้ประมวลผลหรือดำเนินการทางเมตริกซ์ต่อ ไม่ว่าจะนำไป บวก ลบ หรือคูณเมตริกซ์ หรือแม้แต่จะนำไปหาคำตอบของระบบสมการดังที่เคยนำเสนอไปแล้ว ในบทความต่อไปจะนำ ตัวอย่างการใช้คอนสตรัคเตอร์ดังกล่าวในการแก้ปัญหาระบบสมการไม่เชิงเส้น โดยใช้ระเบียบวิธี Newton-Raphson 

เขียนโปรแกรมด้วย C++ : ฟังก์ชันสำหรับโหลดไฟล์ input เพื่อสร้างเตริกซ์

การเขียนโปรแกรมในบทความนี้จะนำเสนอ การสร้างฟังก์ชันสำหรับโหลดไฟล์ input เพื่อสร้างเมตริกซ์ ซึ่งจะนำความรู้เกี่ยวกับการเขียนโปรแกรมการอ่านไฟล์ชนิด text file ที่ได้กล่าวถึงไปแล้วมาประยุกต์กับการสร้าง Dynamics Array แบบ 2 มิติ เพื่อสร้างเมตริกซ์ไว้สำหรับใช้งานต่างๆต่อไป ลักษณะการติดต่อกับผู้ใช้จะเป็นการอ่านข้อมูลจาก text file ซึ่งเราได้กำหนดรูปแบบของการบันทึกข้อมูลไว้เรียบร้อยแล้วนั่นเอง โดยเราจะเรียก text file นั้นว่า input file สมมุติให้นามสกุลไฟล์เป็น .txt โดยเราสามารถใช้โปรแกรม Notepad สำหรับการบันทึกข้อมูล เรามากำหนดรูปแบบของ input file กันก่อนครับ 

ในบรรทัดแรก จะเป็นการบันทึกจำนวนแถวและหลักของเมตริกซ์
ในบรรทัดถัดไปจะเป็นสมาชิกในเมตริกซ์ในตำแหน่งต่างๆ ดัง ตย
3,4
1, 3, 5, 6,
3, 4, 5, 7,
5, 6, 7, 8,

การทำงานของโปรแกรมจะเริ่มจาก
1. รับตำแหน่งและชื่อ input file ที่จะเข้าไปอ่าน
2. อ่านข้อมูลขนาดของเมตริกซ์ แถวxหลัก
3. สร้าง Dynamic Array ตามขนาดของเมตริกซ์
4. อ่านข้อมูล Aij จาก input file ครั้งละ 1 ตัว จนหมด
5. พิมพ์เมตริกซ์ที่อ่านได้จาก input file บนหน้าจอ

ฟังก์ชันอ่านข้อมูลจาก input file มีรูปแบบดังนี้
ReadInput(จำนวนแถว , จำนวนหลัก ,input file name) โดยจะคืนค่าเป็น Dynamics Array แบบ 2 มิติ

ฟังก์ชั่น printMatrix

เรามาดู main Program สำหรับเรียกใช้ฟังก์ชันทั้งสองกันครับ

เมื่อคอมไพล์โปรแกรมนี้จะปรากฎหน้าต่างโปรแกรม

ในบรรทัด Enter path of input file : ให้ท่าน input ชื่อ input file ที่ต้องการอ่านลงไป ตัวอย่างเช่น c:\\aa.txt แล้ว Enter หากไฟล์นั้นมีอยู่จริงและมีข้อมูลที่ถูกต้องโปรแกรมจะแสดงขนาดของเมริกซ์และสมาชิกของเมตริกซ์ให้ทราบดังแสดงในภาพ

ภาพหลังการอ่านข้อมูลจาก text file 

จากตัวอย่างการอ่านข้อมูลจาก text file ดังกล่าว หากเราต้องการดำเนินการระหว่างเมตริกซ์เช่น การบวก ลบ การคูณ การ transpose และให้บันทึกผลการดำเนินการลง text file เพื่อนำไปใช้งานต่อก็สามารถทำได้โดยง่ายแล้ว เดี๋ยวมาเขียนโปรแกรมกันต่อในบทความต่อไปครับ

เขียนโปรแกรมภาษา c++ : ตัวอย่างการบวกและคูณ เมตริกซ์

เขียนโปรแกรมภาษา c++ ในบทความนี้จะขอยกตัวอย่างการประยุกต์ใช้ Dynamics Array และการเรียกใช้ฟังก์ชันมาช่วยในการดำเนินการบวกและคูณเมตริกซ์ ซึ่งเป็นตัวอย่างแสดงการเรียกใช้ฟังก์ชั่นต่างๆ ก่อนที่บทความต่อไปจะนำเสนอ การเขียนโปรแกรมภาษา c++ สำหรับดำเนินการทางเมตริกซ์แบบสมบูรณ์ สามารถนำไปใช้งานได้โดยง่ายครับ
ในลำดับแรกของการเขียนโปรแกรมนี้ เราจะสร้างตัวแปรโครงสร้างขึ้นมาชนิดหนึ่งประกาศไว้ตอนต้นของโปรแกรม นั่นคือตัวแปรโครงสร้าง Matrix ซึ่งภายในประกอบด้วยตัวแปรที่เก็บค่าจำนวนแถว จำนวนหลักและตัวแปรพอยเตอร์สำหรับเก็บสมาชิกในเมตริกซ์ ดังนี้ครับ
struct Matrix
{
int r ; // rows
int c; // columns
float **m;  // data
};

ลำดับถัดมาเป็นการสร้างฟังก์ชันสำหรับใช้ในการดำเนินการระหว่างเมตริกซ์ โดยในตัวอย่างนี้จะมีทั้งหมด 5 ฟังก์ชั่นดังนี้
1. ฟังก์ชั่น creatematrix ทำหน้าที่สร้าง matrix ขึ้นมา โดยจะคืนค่าเป็นตัวแปรโครงสร้างชนิด Matrix ซึ่งจะรับค่าจำนวนแถว(iRow) และ จำนวนหลัก(iCol) ของ Matrix ที่ต้องการสร้างรูปแบบดังนี้
Matrix creatematrix(int iRow,int iCol)

2. ฟังก์ชั่น plusmatrix ทำหน้าที่ดำเนินการบวก matrix สองเมตริกซ์ โดยจะคืนค่าเป็นตัวแปรโครงสร้างชนิด Matrix ซึ่งจะรับค่าเป็น matrix 2 matrix ที่ต้องการนำมาดำเนินการบวกกัน มีรูปแบบดังนี้
Matrix plusmatrix(Matrix a , Matrix b)

3. ฟังก์ชั่น mulMatrix ทำหน้าที่คูณเมตริกซ์ด้วยค่าคงที่ โดยจะคืนค่าเป็นตัวแปรโครงสร้างชนิด Matrix ซึ่งจะรับค่า 2 ค่าคือ ค่าคงที่ที่นำมาคูณ และ matrix ที่ต้องการมาดำเนินการคูณด้วยค่าคงที่ มีรูปแบบดังนี้
Matrix mulMatrix(float fm ,Matrix a)

4. ฟังก์ชั่น Product ทำหน้าที่ดำเนินการคูณระหว่างเมตริกซ์ โดยจะคืนค่าเป็นตัวแปรโครงสร้างชนิด Matrix ซึ่งจะรับค่าเมตริก 2 เมตริกซ์ ที่ต้องการมาดำเนินการคูณระหว่างเมตริกซ์ มีรูปแบบดังนี้
Matrix Product(Matrix a , Matrix b)

5. ฟังก์ชั่น freedynamicarray ทำหน้าที่คืนหน่วยความจำให้ระบบหลังการใช้งาน Dynamic Array ทั้งหมดของโปรแกรม

เรามาดู Main Program ของการเขียนโปรแกรม แสดงการบวกและคูณเมตริกซ์กันครับ

มาดูผลลัพธ์จากการคอมไพล์ครับ

ผลการคอมไฟล์โปรแกรม การคูณกันระหว่างเมตริกซ์

จากตัวอย่างการเขียนโปรแกรมที่ได้นำเสนอมา จะพบว่าการที่จะเขียนโปรแกรมให้สมบูรณ์จะต้องเพิ่มโค้ดในการตรวจจับความผิดพลาดของการคำนวณทางตัวเลขเพื่อป้องกันความผิดพลาดที่จะเกิดขึ้น จากตัวอย่างนี้เราจะพัฒนาโปรแกรมให้มีความสมบูรณ์มากขึ้นโดยจะเพิ่มการลบเมตริกซ์ การสร้างเมตริกซ์มาตรฐานต่างๆ ตัวดำเนินการแยกองค์ประกอบเมตริกซ์ และสุดท้ายการกำหนดสมาชิกในเมตริกจากการอ่านข้อมูลรูปแบบ text file ที่ผมได้นำเสนอไปแล้วในบทความก่อนหน้านี้ ซึ่งจะทำให้การใช้งานโปรแกรมนี้สามารถใช้งานได้ง่ายที่สุด นอกจากนี้จะทำการพัฒนา Class CMatrix เพื่อให้ผู้ที่สนใจนำไปประยุกต์ใช้งานกับการสร้าง Application ต่างๆ ต่อไป โดยจะเพิ่มความสามารถในการแก้ระบบสมการเพื่อนำไปใช้ในการแก้ปัญหาทางวิศวกรรมและระเบียบวิธีไฟไนต์เอลิเมนต์ต่อไปครับ

หัวข้อที่น่าสนใจ

เขียนโปรแกรม c++ เพื่อสร้างฟังก์ชันสำหรับโหลด Input File และสร้างเมตริกซ์ 

การประยุกต์ใช้ Class CMatrix ในการสร้างโปรแกรมช่วยแก้ระบบสมการเชิงเส้น