Matlab with LaTeX interpreter

เคยต้องการพล็อตกราฟของสมการคณิตศาสตร์ด้วย Matlab ไหมครับ สิ่งหนึ่งที่น่าปวดหัวคือการผนวกสมการลงไปในกราฟด้วย เช่นอาจจะต้องการพล็อตกราฟ $$e^{2\log x}$$ แล้วต้องการเขียนสมการกำกับที่กราฟ แบบสิ้นคิดทำได้ดังนี้
จะเห็นได้ว่ารูปข้างบนไม่ต้องทำอะไรมาก ใช้คำสั่ง

 text(2,2,'exp(2log(x))')

เป็นอันเสร็จพิธี

แต่การทำอะไรให้มันเต็มความสามารถนั้น ถึงแม้จะไม่ได้คะแนนเพิ่ม แต่ทำให้เอกสารอ่านง่ายดูดี เป็นเรื่องที่ควรกระทำ ดังนั้นใช้ความรูป LaTeX เพิ่มเข้าไปดังนี้

text(2,2,'$e^{2\log(x)}$','interpreter','latex','fontsize',18)

ซึ่งจะได้ผลของสมการที่สวยงามดังรูป
ซึ่งวิธีการนี้ใช้ได้หมดนะครับ ไม่ว่าจะเป็นคำสั่ง xlabel, ylabel, และ title แต่ช้าก่อนมันใช้กับคำสั่ง legend ไม่ได้ และนี่เป็นที่มาของบทความนี้ ในกรณีที่เราต้องการใช้คำสั่ง LaTeX กับ legend โดยสั่งจาก Matlab โดยตรก ไม่ต้องไปใช้เมาส์คลิ๊กขวาแล้วเลือก เราจะทำอย่างไร สมมติว่าผมพล็อตกราฟ $$x^3$$ และ $$x^2$$ บนแกนเดียวกัน แล้วใช้คำสั่ง

legend('x^3','x^2');

ผมก็จะได้กราฟตามรูป

ซึ่งก็จะได้ $$x^3$$ และ $$x^2$$ ที่ดูดีพอสมควร (กรณีนี้คำสั่ง LaTeX บางคำสั่ง รวมไปถึงอักษรกรีกบางตัวจะใช้ไม่ได้) แต่ถ้าต้องการให้ $$x^3$$ และ $$x^2$$ สวยงามขึ้น เราไม่สามารถใช้วิธีการตามข้างบนได้ จะต้องทำดังนี้

s =  legend('$x^3$','$x^2$');
set(s,'interpreter','latex');

ซึ่งก็จะได้กราฟสวย ๆ ตามต้องการครับ ดังแสดงในรูปข้างล่าง

หวังว่าจะมีประโยชน์ไม่มากก็น้อยนะครับ อนึ่งรูปที่โพสลงไปเป็น png คุณภาพอาจจะไม่ดีนัก แต่ของจริงดูดี รับรอง

Advertisements

การดึงข้อมูลจากกราฟใน Matlab

ใช้ Matlab แล้วเคยประสบปัญหาแบบนี้ไหมครับ

“ต้องการพล็อตกราฟใหม่ให้เหมือนรูปเก่า แต่หา m-file ไม่เจอ”

และมีปัญหาต่อเนื่องมา

“มีรูปสองรูปที่พล็อตไว้แล้วต้องการเปรียบเทียบบนแกนเดียวกัน แต่ข้อมูลที่พล็อตไว้ค่าบนแกน x ไม่เหมือนกัน”

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

โม้ซะเยอะ ถ้าเรามีกราฟ Matlab ในมือ แล้วทำ m-file ที่ใช้งานหาย หรือหาไม่เจอ แล้วต้องการนำข้อมูลจากกราฟที่พล็อตไปพล็อตเปรียบเทียบกับข้อมูลอื่น ในคู่มือ Matlab นั้นมีวิธีการ(ซ่อน)อยู่ในคู่มือ ซึ่งน้อยคนนักที่จะหาเจอ ผมเจอวิธีการนี้ตั้งแต่ปี 2007 ได้จดบันทึกเอาไว้ วันนี้เอามาเขียนใหม่แล้วกันครับ

ให้เลือกกราฟที่เราต้องการ จากนั้นให้ใช้คำสั่ง findobj เพื่อสร้างตัวแปรวัตถุของกราฟนั้น ๆ จากนั้นก็ดึงข้อมูลในแกน x และ y มาใช้ได้เลย ตามคำสั่งข้างล่าง

pts=findobj(gca,'type','line');
x = get(pts,'XData');
y = get(pts,'YData');

ในกรณีมีข้อมูลมากกว่าหนึ่งชุด เราสามารถกำหนดเลขดัชนีสำหรับข้อมูลในกราฟได้ดังนี้

y1 = get(pts(1),'YData');
y2 = get(pts(2),'YData');

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

Hold all

ทุกคนที่เคยใช้ Matlab คงจะไม่มีใครไม่เคยใช้ฟังก์ชัน hold on หรือ hold off  ที่ blog ของ Loren มีการนำเสนอตัวเลือกของฟังก์ชัน hold นั่นคือ hold all และเหมือนทุกครั้ง ผมก็พึ่งรู้ว่ามีงี้ด้วย

ยกตัวอย่างปัญหาที่มักเกิดขึ้นเสมอ ๆ แล้วกันครับ เวลาเราพล็อตกราฟเพื่อเปรียบเทียบ เอาตามตัวอย่างของ Loren แต่นำเสนอให้เข้ากับสถานการณ์จริงมากขึ้นคือ สมมติว่า  เราสร้างฟังก์ชัน sin แล้วพล็อตดังนี้

t = 0:0.005:1;
f = sin(2*pi*10*t);
plot(t,f)

ซึ่งเราก็จะได้กราฟดังรูปsin

Continue reading

จงเชื่อคำเตือนของ MATLAB

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

ผมจำเป็นต้องวาดกราฟของการลู่เข้าของค่าพารามิเตอร์ซึ่งมีประมาณ 2 หมื่นค่า อาจจะสร้างเล่น ๆ ได้ดังนี้

theta = [];

for ii = 1:20000
theta(1:20,ii) = rand(20,1);
end

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

theta = [];

theta = zeros(20,20000);
for ii = 1:20000
theta(1:20,ii) = rand(20,1);
end

ลองรันชุดคำสั่งที่สองดู ยังไม่ทันกระพริบตาก็เสร็จแล้ว ซึ่งชุดคำส่งแรกใช้เวลาประมาณ 170 วินาที เกือบ ๆ สามนาที ในขณะที่ชุดคำสั่งที่สองใช้เวลาเกือบ 0.1 วินาที ซึ่งต่างกันมากทีเดียว

ก็ไม่แปลก เพราะเวลาไม่ได้กำหนดขนาดของอะเรย์ เวลามีการขยายหรือทำอะไรกับมัน ย่อมน่าจะมีกระบวนการต่าง ๆ พอสมควร (ผมไม่รู้และไม่สนใจ) แต่ที่แปลกคือมันต่างกันขนาดนี้เลยเหรอ ดังนั้นแล้วเวลาใช้งาน MATLAB ก็ควรจะใส่ใจต่อคำเตือนของ IDE ของ MATLAB ด้วย

Why shouldn’t use determinant to check the singularity of matrices

เป็นที่รู้กันมาตั้งแต่ชั้นมัธยมแล้วว่าเราจะหาส่วนกลับของเมตริกซ์ได้ก็ต่อเมื่อค่า determinant ไม่
เท่ากับเมตริกซ์ 0 พอเรียนสูงขึ้นเขาก็บอกว่าไม่ให้ทำแบบนี้นะ เพราะสภาวะทางตัวเลขที่คำนวณ
ด้วย determinant นั้นไม่ดี ให้ตรวจสอบด้วย Rank และการหา Rank ของเมตริกซ์ที่ดีที่สุดก็ให้ใช้
Singular Value Decomposition ซึ่งเมตริกซ์ที่จะหาส่วนกลับได้นั้นจะต้อง full Rank

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

     3825       -1713       -1061        -497         710        -502       -3471      3136
    -1713        3834       -1632         161        -527         722        -412      -3561
     -1061       -1632        4563         983        -109       -419       4004      -1222
      -497         161         983         315        -178          23           832     -978
        710        -527        -109        -178         415        -218         -29       1132
      -502         722        -419          23        -218         431          55        -149
     -3471        -412        4004         832         -29          55        6640     -1265
      3136       -3561       -1222        -978        1132        -149     -1265     7539

ค่า determinant ของตัวเลขชุดนี้คือ –54413 นั่นหมายความว่าเราต้องหาส่วนกลับของเมตริกซ์
ตัวนี้ได้สิ แต่ไม่ใช่เลย ถ้าใช้คำสั่ง rank ใน MATLAB เราจะได้ 7 ซึ่งเมตริกซ์ไม่ full rank ถ้าใช้
svd หา จะได้ rank เป็น 6 หมายความว่าเมตริกซ์นี้หาส่วนกลับไม่ได้

เชื่อแล้วครับผม

ปล. ถ้าคิดด้วยมือ determinant จะไม่ผิดนะครับ

function reshape ของ MATLAB

เป็นฟังก์ชันที่ถ้าใช้บ่อยแล้วจะดี เพราะการกระทำใด ๆ ในรูปเมตริกซ์โดยการหลีกเลี่ยงการใช้
การวนลูปแล้วกระทำกับตัวเลขทีละตัวจะทำให้การคำนวณโดยใช้ MATLAB เร็วขึ้น

ปัญหาของผมคือผมมีข้อมูล

reshape1

ผมต้องการให้เป็นแบบนี้

reshape2 ในคู่มือของ MATLAB บอกว่าให้ใช้คำสั่ง reshape ซึ่งมีวิธีการใช้ดังนี้

B = reshape(A,m,n)
โดยที่ A คือเมตริกซ์ที่ต้องการจะทำการเปลี่ยนแปลง
m และ n คือจำนวนแถวและหลัก

ถ้าผมใช้ B = reshape(y,[],1)  ผมจะได้

reshape3.

ซึ่งไม่ใช่สิ่งที่ต้องการ คิดอยู่ต้องนานว่าจะทำอย่างไร สรุปว่าก็แค่ B = reshape(y’,[],1) นั่นก็
คือทำการ transpose ก่อนนั่นเอง

 

 

วาดวงรีด้วย MATLAB

จากความตอนที่แล้ว

ขั้นต่อไปคงเป็นการสร้างวงรี จริง ๆ แล้วเราเรียนเรื่องวงรีจากเรขาคณิตวิเคราะห์ในบทที่ว่าด้วย
ภาคตัดกรวยที่มี วงกลม วงรี พาราโบลา ไฮเปอร์โบลา ยากสุด ๆ กลับไปดูในแบบเรียนเอาเอง
นะครับ หรือที่ Conic Section Gallery

การวาดวงรีด้วย MATLAB สำหรับผมแล้วผมจะไม่ใช้เรขาคณิตวิเคราะห์ เพราะผมเป็นวิศวกรที่
ไม่เขียนแบบด้วยไม้บรรทัดและวงเวียนอีกแล้ว ผมจะวาดแบบง่ายโดยให้ค่าบนแกน y = r_1\sin\theta และค่าบนแกน x = r_2\cos\theta ซึ่งคราวนี้ค่า r_1 และ r_2
จะมีขนาดไม่เท่ากัน

รหัสต้นของ MATLAB ที่ดูดีขึ้นก็จะเป็น

theta = linspace(0,2*pi,100);
r = 1;
x = 2*r*cos(theta);
y = r*sin(theta);
axis equal
axis([-3,3,-3,3]);

el1

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

g= U_\theta u, \qquad U_\theta = \begin{bmatrix}\cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix}

โดยที่ u = [x \quad y]^T คือค่าบนแกนนอนและแกนตั้งตามลำดับ ก็จับมาคูณ
กันดื้อ ๆ เลยครับ ได้รหัสต้นดังนี้

ar = 30*pi/180;
Uth = [cos(ar) -sin(ar); sin(ar)  cos(ar)];
g = Uth*[x ; y ];
plot(g(1,:),g(2,:));
axis([-3,3,-3,3]);

el2

ก็จะได้รูปวงรีที่เอียงทำมุม 30 องศากับแกนนอนตามต้องการ

อนึ่ง เพื่อให้ได้อะไรบ้าง U_\theta จะทำการเปลี่ยน orthogonal basis
จาก e_1 = [1 \quad 0]^T และ e_2 = [0\quad 1]^T ไปเป็น e_1 = [\cos\theta \quad \sin\theta]^T
และ e_2 = [-\sin\theta \quad \cos\theta]^T นั่นเอง

การสร้างวงกลมด้วย MATLAB

เห็นมีคำค้นการสร้างวงกลมด้วย MATLAB ไม่ทราบว่าคนค้นได้คำตอบหรือยัง ก็ถือ
โอกาสเขียนวิธีการไว้ตรงนี้เลย เผื่อว่าคนที่ค้นจะกลับมาอีกแล้วจะได้คำตอบ

ว่ากันตามตรงการสร้างวงกลมด้วย MATLAB ถือเป็นปัญหาเส้นผมบังภูเขานะครับ
สมการวงกลมนั้นมีลักษณะดังนี้

x^2+y^2 = r^2

โดยที่ r คือรัศมีของวงกลม ถ้าผู้ต้องการวาดวงกลมด้วย MATLAB เอาเส้นผมออก
ก็จะพบว่าจริง ๆ แล้วเราสามารถเขียนสมการข้างต้นได้ดังนี้

r\sin^2\theta + r\cos^2\theta = r

โดยให้ค่าบนแกน y คือ r\sin\theta  และค่าบนแกน x คือ r\cos\theta
แค่นี้เราก็สามารถนำมาสร้างคำสั่งของ MATLAB ได้แล้วง่าย ๆ ดังนี้

theta = 0:0.1:2*pi;
r = 1;
y = r*sin(theta);
x = r*cos(theta);
plot(x,y);

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

MATLAB Central – Loren on the Art of MATLAB

บทความอาทิตย์ที่แล้วของเธอกล่าวถึงเรื่องของฟังก์ชัน hypot ว่ามีไว้ทำไม ก็พึ่ง
รู้จักฟังก์ชันนี้เหมือนกัน

HYPOT Robust computation of the square root of the sum of squares
C = HYPOT(A,B) returns SQRT(ABS(A).^2+ABS(B).^2) carefully computed to avoid underflow and overflow.

ข้างบนเป็น Help ที่ได้จากตัว MATLAB

เธอบอกว่าสำหรับตัวเลขทั่วไปแล้วไม่มีปัญหาในการคำนวณค่าของฟังก์ชัน

y = \sqrt{a^Ta + b^Tb}

นี้หรอก

ที่น่าสนใจคือฟังก์ชันนี้จะให้ผลลัพท์เป็น Inf เมื่อตัวเองเกินค่า realmax ซึ่งถ้า
เป็นวิธีการอื่น ๆ จะได้ค่า realmax แทนที่จะเป็น Inf

บทความของเธอมีหลายอันน่าสนใจดีนะครับ สำหรับคนชอบ MATLAB

เขียน MATLAB ให้เร็วขึ้น

สำหรับ MATLAB นั้น มันเป็นกรอบการทำงานที่ครอบจักรวาลเกินไป ตัวใหญ่และ
อ้วน แม้แต่คนของ Mathworksเองก็รู้อะไรไม่เยอะ การเขียน m-files ให้เร็วขึ้น
จึงเป็นเรื่องพึงปรารถณา (Toolbox ควรจะเขียนเป็น m-files เพราะพัฒนาต่อยอด
ได้ง่าย การทำงานแบบวิเคราะห์หรือสังเคราะห์ ความเร็วไม่ใช่เรื่องจำเป็นนัก Prof.
S. Boyd เนี่ยแหละบอกว่ามันรันช้านัก ก็ไปเที่ยวปีหนึ่งแล้วค่อยกลับมาทำงานต่อ
มันก็เร็วแล้วหล่ะ) เพราะการทำงานในขั้นตอนพัฒนานั้นต้องการความเร็ว เพื่อให้
ผู้พัฒนาสามารถพัฒนาได้เร็วขึ้น ขืนรัน MATLAB เพื่อแสดงผลทีข้ามวันย่อมไม่
ใช่ผลดี วันนี้มีตัวอย่างการเขียน MATLAB ให้เร็วขึ้น ดูตัวอย่างนี้ครับ

D(t) =\displaystyle{\left\{\sum^{n_a}_{i=1}[a_i-\hat{a}_i(t)]^2+ \sum^{n_b}_{i=0}[b_i-\hat{b}_i(t)]^2\right\}^{1/2}}

ตัวอย่างนี้ดีมาก ๆ

D = sqrt((sum(a-ah(t)))^2 + (sum(b-bh(t)))^2);

ง่าย ๆ แบบนี้

ถ้าใช้ความรู้ทางพีชคณิตเชิงเส้นก็จะทำแบบนี้

D = sqrt((a-ah(t))’*(a-ah) + (b-bh(t))’*(b-bh(t)));

วิธีไหนที่เร็วกว่ากัน คิดกันง่าย ๆ ถ้า a, b เป็นเวกเตอร์ที่มีสมาชิก n ตัว แบบที่
หนึ่งจะมีการบวกกัน 2n +1 ครั้ง และมีการยกกำลังสองหรือคูณ 2 ครั้ง นั้นคือมีการกระทำ
ทางพีชคณิต 2n+3 ครั้งแล้วค่อยถอดรากที่สอง ในขณะที่แบบที่สอง มีการคูณกัน 2n ครั้ง มี
การบวกกัน 1 ครั้ง นั้นคือมีการกระทำทางพีชคณิต 2n+1 ครั้ง (ไม่นับการ transpose ซึ่ง
เสียเวลาเหมือนกัน) ไม่น่าเชื่อแบบที่สองช้ากว่าแบบแรกเกือบสองเท่าทีเดียว (เมื่อปริมาณ
ข้อมูลเพิ่มขึ้น)