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