jump to navigation

Numerical methods: Solving ODE October 8, 2010

Posted by viboon in : Mathematics , add a comment

Given

%Initial conditions
x=0;
y=1;
%Step size
h=1; %Except Runge-Kutta methods
%Span
n=20;

Euler’s method

for i=1:1:n
   y = y + f(x,y)*h;
   x = x + h;
end

Modified Euler’s method

for i=1:1:n
    y1 = y + f(x,y)*h/2;
    x1 = x + h/2;
    y = y + f(x1,y1)*h;
    x = x + h;
end

Heun’s method

for i=1:1:n
   y1 = y + f(x,y)*h;
   x = x + h;
   a = (f(x,y) + f(x,y1))/2;
   y = y + a*h;
end

3rd-order Runge-Kutta method

h = 0.25;
for i=1:1:n
   f1 = f(x,y);
   x0 = x + h/2;
   y0 = y + f1*h/2;

   f2 = f(x0,y0);
   x0 = x + h;
   y0 = y - f1*h + 2*f2*h;

   f3 = f(x0,y0);
   y = y + (f1 + 4*f2 + f3)*h/6;
   x = x + h;
end

4th-order Runge-Kutta method

h = 0.25;
for i=1:1:n
   f1 = f(x,y);
   x0 = x + h/2;
   y0 = y + f1*h/2;

   f2 = f(x0,y0);
   y0 = y + f2*h/2;

   f3 = f(x0,y0);
   x0 = x + h;
   y0 = y + f3*h;

   f4 = f(x0,y0);
   y = y + (f1 + 2*f2 + 2*f3 + f4)*h/6;
   x = x + h;
end

Tags: , , , , ,

Related posts:

ODE เบื้องต้นสำหรับ MATLAB August 19, 2010

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

Numerical methods นิยมใช้ในการแก้ปัญหาแบบ initial-value problem โดยอาศัยการประมาณเพื่อให้ได้ค่าที่ใกล้เคียงกับ exact solution

ใน MATLAB สามารถแก้ปัญหา ODEs ได้โดยใช้ functions ทาง ODE ต่างๆ ที่มีให้ ซึ่งจะแตกต่างไปตามความแม่นยำและระยะเวลาในการ solve รูปแบบการใช้งานเป็นไปตามคำสั่งด้านล่าง

1
2
[outputs] = dydt(inputs) %Derivative function(s)
[t,y] = solver(@dydt,TimeInterval,InitialConditions,options)

โดย function ที่นิยมใช้ในการคำนวณ (solver) มีดังนี้
ode45 — ความถูกต้อง ปานกลาง (โดยทั่วไปจะใช้ solver ตัวนี้)
ode23 — ความถูกต้อง ต่ำ (น้อยกว่า ode45)
ode113 — ความถูกต้อง สูงไปจนถึงต่ำ (เหมาะกับปัญหาที่มีความซับซ้อน)
ode15s — ความถูกต้อง ปานกลางไปจนถึงต่ำ (มักใช้เมื่อ ode45 ไม่สามารถทำการคำนวณได้ถูกต้องหรือให้ผลการคำนวณที่แกว่งมาก (stiff problems))

สำหรับการเลือกใช้วิธีในการแก้ปัญหาก็แตกต่างไปตามความถูกต้องและความเร็วในการคำนวณ โดยจะขึ้นอยู่กับ order ของ derivative function เช่น
Euler’s Method – 1st order expansion
Midpoint method – 2nd order expansion
Runge-Kutta – 4th order expansion


Tags: , ,

Related posts:

แก้ปัญหา 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: