jump to navigation

แก้ปัญหา Ordinary Differential Equation ด้วย MATLAB August 13, 2010

Posted by viboon in : Mathematics, Matlab , add a comment

ยกตัวอย่างการแก้ปัญหาที่เป็น second-order system อย่างระบบสปริงที่มีการเคลื่อนที่แบบ 1 มิติในแนวราบ (ไม่คิดแรงโน้มถ่วง)

F=m*x”(t)+k*x(t)

โดยมี Initial conditions คือ
x(0)=0 และ x'(0)=1
กำหนดให้ m=1 และ k=5 หากไม่คิดแรงภายนอก F=0

เนื่องจากสมการนี้สามารถปรับให้อยู่ในรูปของ first-order ได้ โดยกำหนดให้
y=x’ และ y’=x”
ดังนั้น สามารถเขียนสมการได้ใหม่ว่า

0=m*y'(t)+k*x(t)

จากตรงนี้เราสังเกตได้ว่ามีสมการอยู่ 2 ชุด คือ spring equation และ transform equation ซึ่งก็คือ

0=m*y'(t)+k*x(t) และ
y=x’

จัดรูปใหม่ให้เป็นรูปแบบเดียวกัน เป็น

y’=-k*x(t)/m
x’=y

เพื่อความสะดวกในการคำนวณ เราสามารถเขียน 2 สมการนี้ให้อยู่ในรูปของ vector หรือ matrix คือ z=[x;y] หรือก็คือ vector ที่บอกตำแหน่ง (x) และความเร็ว (x’) นั่นเอง เมื่อเขียน z ให้อยู่ในรูป z’ ก็จะได้ z’=[x’;x”] หรือก็คือ z’=[y;y’]

ดังนั้น
z’=[z(2);-k*x(t)/m] โดยที่ z(2) ก็คือ y นั่นเอง

ทีนี้เราก็สามารถ solve ปัญหา second-order ในรูปของ first-order ได้โดยการใช้ transform equation ซึ่งก็คือ y=x’ ได้

1
2
3
4
5
6
7
8
m=1;
k=5;
initcon=[0; 1]; %Initial conditions
%กำหนดให้ spring เป็น function ของ z'
%ซึ่ง z'(1)=z(2)=y=x' และ z'(2)=y'=x''
spring = @(t,z)[z(2); -k * z(1)/ m]; %Using anonymous function
tspan=[0 10]; %Time span
[t,z] = ode23(spring, tspan, initcon);

สังเกตได้ว่าค่าที่ได้คือ เวลา และ vector z ซึ่งในที่นี้ z จะเป็นผลลัพธ์ของการคำนวณ ซึ่งก็คือ z=[x;y] หรือ z=[x;x’] นั่นเอง หรือก็คือ ตำแหน่ง และ ความเร็วของมวล

1
2
3
4
subplot(2,1,1);
plot(t,z(:,1)), title('Position vs. Time');
subplot(2,1,2);
plot(t,z(:,2)), title('Velocity vs. Time');

ดูรายละเอียดการใช้ ODE ด้วย MATLAB ได้ที่:
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/ode23tb.html


Tags: , , ,

Related posts: