Phase Portraits with MATLAB


ทุกครั้งที่เริ่มเรียน Nonlinear System สิ่งที่เรียนต้องรู้ก็คือการตรวจสอบเสถียรภาพ
ของระบบ และสิ่งที่จะถูกหยิบยกมาเป็นตัวอย่างเสมอ ๆ ก็คือระบบลำดับที่สอง จาก
นั้นก็ใช้ Phase Portraits วิเคราะห์ระบบ ก็จะมีรูปโชว์ แล้วก็สอนวิธีการวาดด้วยมือซึ่ง
ใช้ได้กับระบบง่าย ๆ เท่านั้น สรุปง่าย ๆ ก็คือเรียนวิธีวาดไว้เพื่อออกข้อสอบ พอไปใช้
งานจริงก็ต้องพึ่ง MATLAB และความจริงอันโหดร้ายก็คือไม่มีใครสอนวิธีสร้างรูปพวก
นี้ในชั้นเรียนวิชา Nonlinear Systems ด้วย MATLAB ส่วนใหญ่ก็จะเอารูปจากใน
หนังสือมาสอน เมื่อก่อนราว ๆ 7 ปีที่แล้ว ผมได้เรียนวิชานี้ผมต้องพึ่งพาอาศัยหนังสือ
เล่มนี้ (ซื้อเก็บไว้ตอนเป็นเวอร์ชัน 2) Mastering SIMULINK ซึ่งเปิดเผยทุกอย่าง
ปัจจุบันไม่ต้องแล้วครับ หนังสือแนว How To เหล่านี้นอกจากจะล้าสมัยแล้วยังราคา
แพง ทางที่ดีก็หาคำค้นแล้วอาศัย Google เอาดีกว่า แต่บางทีการหาจาก Google ก็
ใช่ว่าดี ทางที่ดีจดใส่บล็อกดีกว่าครับ อย่างเรื่องนี้เป็นต้น

สำหรับใครที่ต้องการใช้งานโปรแกรมที่สามารถพล็อตสนามเวกเตอร์หรือ Phase
plane ได้ ก็ไปโหลดโปรแกรมที่นี่ได้เลย ODE software for MATLAB มีทุกอย่าง
ให้เสร็จ แต่ถ้าใครต้องการรายละเอียดที่อธิบายเข้าใจก็ที่นี่ Teaching material of
Frank Schilder
ซึ่งผมจะสรุปเอกสารของแกไว้ย่อ ๆ ดังนี้

กรณีที่ต้องการ Phase Portraits ของ Nonlinear systems ที่อธิบายด้วยสมการข้าง
ล่าง ซึ่งเป็นสมการเชิงพลวัตรของการแกว่งของลูกตุ้ม
\dot{x}_1 = x_2; \; \; \dot{x}_2 = -\theta^2\sin x_1 - bx_2
สำหรับระบบง่าย ๆ การสร้าง anonymous ดูจะเหมาะสมกว่าดังนี้

b = 0; theta = 1;
f = @(t,x) [ x(2); -theta^2*sin(x(1)) – b*x(2)

เวลาจะแก้สมการอนุพันธ์ก็ใช้คำสั่ง

ode45(f, [0 100], [0.1 0])

พารามิเตอร์ที่ใส่ก็น่าจะชัดเจน ถ้าอยากเปลี่ยนแปลงแก้ไขอะไรก็อย่าลืม help ode45

ในการพล็อต phase portrait สิ่งที่สำคัญอีกอันหนึ่งคือ direction field ซึ่งทำได้
โดยใช้ฟังก์ชัน quiver โดยทำได้ดังนี้ (เอาตัวอย่างมาจากเอกสารของ Schilder)

[xx yy] = meshgrid(-1.5*pi:pi/10:1.5*pi, -pi:pi/10:pi);
fx = []; fy = [];
for i=1:size(xx,1)
for j=1:size(xx,2)
ff = f(0,[xx(i,j) yy(i,j)]);
fx(i,j)=ff(1); fy(i,j)=ff(2);
end
end
quiver(xx, yy, fx, fy);
axis([-1.5*pi 1.5*pi -pi pi]);

ซึ่งจะได้รูป direction field (vector filed) ตามรูปข้างล่าง จริง ๆ การพล็อตรูปนี้ถ้า
อ่านคู่มือของคำสั่ง quiver ก็สามารถทำได้อยู่แล้ว

ในการวาด phase portrait นั้นจะต้องใช้ความรู้เรื่อง ODE ซึ่งต้องแก้สมการเชิง
อนุพันธ์ โดยกำหนดค่าตั้งต้นที่เหมาะสม โดยแก้สมการตามเส้นแนวดิ่ง x=0 (เกี่ยวกับ equilibrium points) โดยการวาดครึ่งแรกด้วยรหัสดังนี้

y0 = -pi:pi/10:pi;
for i=1:length(y0)
[t,z] = ode45(f, [0 10], [0 y0(i)]);
hold on;
plot(z(:,1), z(:,2), ’r-’, ’LineWidth’, 2.0);
drawnow;
hold off;
end

อีกครึ่งก็เปลี่ยนฟังก์ชันเป็น fb = -f

fb = @(t,x) -f(t,x);
y0 = -pi:pi/10:pi;
for i=1:length(y0)
[t,z] = ode45(fb, [0 10], [0 y0(i)]);
hold on;
plot(z(:,1), z(:,2), ‘r-‘, ‘LineWidth’, 2.0);
drawnow;
hold off;
end

จะได้รูปดังนี้

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

Advertisements

9 comments on “Phase Portraits with MATLAB

  1. สวัสดีครับ

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

    ที

  2. ค้นข้อมูลในกูเกิลอยู่เลยมาเจอบล็อกนี้เข้า
    เลยได้เทคนิค การพล็อต phase-plan โดยใช้matlab ให้สวยงาม ได้มารู้จักคำสั่งquiver วาดได้สวยขึ้นเยอะเลยครับ

  3. ดีใจมาก ครั บมากครับที่มาเจอบล็อกนี้
    ผมเองก็ ฝันว่า สักวันหนึ่งจะเป็นเจ้าของบล็อกดีๆ
    แบบนี้บ้างครับ แต่ตอนนี้ยังความรู้น้อย ขออาศัยอ่านไปก่อน
    นะครับ บอกได้คำเดียวว่า ปลื้มสุดๆครับ

  4. ดีใจขนาดนั้น ขอบคุณนะครับที่ปลื้ม

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

    นั่นหมายความว่าผู้อ่านสำคัญกว่าผู้เขียน ผู้อ่านจะเป็นตัวผลักดันให้ผู้เขียน
    ต้องตอบคำถาม แล้วได้ใช้ประโยชน์ร่วมกันนั่นเอง

    บล็อกมันจึงเป็นลักษณะที่เก็บเนื้อหามากกว่า ว่าไปแล้วผมเองใช้
    ประโยชน์จากเว็บผมเองบ่อยมาก ๆ

  5. ขอถามเรื่องการ ใช้คำสั่ง SVD ใน MATLAB น่ะครับ คือ ผมอยากได้ eigenvalue และ eigenvector ซ้าย,ขวา น่ะครับ
    ผมก็ ทำแบบนี้ละครับ
    [A,B,C]=SVD(JR);
    แต่ eigenvalue หรือ (B) ที่ได้มาเรียงจากมากไปน้อย คือเรียงให้เราเลยเสร็จ ทำให้ eigenvector ซ้าย,ขวา ก็จะเรียงสอดคล้องกับ B
    ปัญหามันอยู่ที่ว่า ยังต้องการอ้างอิงตำแห่นงของ JR น่ะครับ(แถงและหลัก) แต่พอ SVD แล้วมาเรียงให้เราใหม่แบบนี้ ผมจะทราบได้อย่างไร ละครับว่า eigenvalue ตัวใหน ตำแห่นงของ JR ครับบบบ
    งง มั๊ย เนี่ยยย กรุณาด้วย ครั๊บบ

  6. ขอบคุณมากครับ
    สรุป คือ ผมใช้ คำสั่ง eig แทนครับ
    [A,B] = eig(JR)
    แล้ว C=inv(A)
    eigenvalue หรือ (B) รวมถึง right eigenvector(A) จากคำสั่ง eig แถวและหลักจะมีความสำพันธ์ กับ JR โดยตรง ไม่มีการสลับแถวและหลักนะครับ ผมคิดว่า จะ จะเป็นวิธี singular transformation นะครับ เป็นคนละวิธีกับ SVD แต่น่า จะใช้แทนกันได้ใช่มะครับ ในกรณีที่ที่ ต้องการ eigenvalue และeigenvector น่ะครับ

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s