โปรแกรมช่วยหาคำตอบระบบสมการไม่เชิงเส้นด้วย 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