jump to navigation

การคำนวณการปริมาตรใต้พื้นผิวด้วย Trapezoidal numerical integration June 28, 2011

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

Trapezoidal numerical integration เป็นการหา Integral ของ function ด้วยการประมาณพื้นที่ใต้กราฟ(หรือพื้นผิว)เป็นส่วนย่อยๆ ที่มีรูปร่างเป็นพื้นที่แบบสี่เหลี่ยมคางหมู (Trapezoid)

ในตัวอย่างนี้กำหนดให้หาปริมาตรใต้ surface function

Z = X^2 \sin \left[ {5\left( {X - Y} \right)} \right]

โดยการหา double integral ผ่านวิธี Trapezoidal method

1
2
3
4
5
6
7
xdata = [0; rand(100,1); 1];
ydata = [0; rand(100,1); 1];
x = sort(xdata);
y = sort(ydata);
[X,Y] = meshgrid(x,y);
Z = X.^2.*sin(5*(X-Y));
trapz(y,trapz(x,Z,2),1)

ในกรณีนี้จะได้คำตอบ

ans =
    0.0497

เมื่อใช้วิธี double integral แบบทั่วไป

1
2
F = @(x,y)(x.^2).*sin(5*(x-y));
dblquad(F,0,1,0,1)

คำตอบคือ

ans =
    0.0498

Tags: , , , , , ,

Related posts:

พื้นฐานการคำนวณทางคณิตศาสตร์ด้วย Python March 17, 2011

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

บวก ลบ คูณ หาร modulus ยกกำลัง ทั่วไป

>>> 2+2
4
>>> 2*3
6
>>> 2**3
8
>>> 4/2
2
>>> 5/2
2
>>> 5.0/2
2.5
>>> 5/2.
2.5
>>> 5%2
1
>>> x=input()
10
>>> 10+x
20
>>> y=5
>>> x=30
>>> y*x
150

การเรียกใช้ functions ทาง math

>>> abs(-7)
7
>>> pow(2,3)
8
>>> import math
>>> math.sqrt(25)
5.0

จะสังเกตว่า การเรียกใช้ function บางตัวต้องทำการ import module ชื่อ math เข้ามาก่อน รูปแบบของการเขียนจะเป็น

import modulename

และการเรียกใช้ function จะใช้คำสังในรูปแบบ

modulename.function()

เช่น

1
2
3
x=25
import math
math.sqrt(x)

หรืออาจเขียนให้สะดวกขึ้นเป็น

1
2
3
4
y=input("Enter y value: ")
import math
myfun=math.sqrt
myfun(y)

Tags: ,

Related posts:

Polyhedra March 11, 2011

Posted by viboon in : Mathematics , add a comment

Euler ค้นพบว่าโครงสร้าง 3D พื้นฐานที่มีด้านแต่ละด้านเหมือนกันสามารถเขียนออกมาเป็นความสัมพันธ์ได้ดังนี้ โดย \chi คือ Euler characteristic และจะมีค่าเท่ากับ 2 เสมอ

 \chi  = V - E + F = 2

Name Vertices
V
Edges
E
Faces
F
Euler characteristic:
VE + F
Tetrahedron 4 6 4 2
Hexahedron or cube 8 12 6 2
Octahedron 6 12 8 2
Dodecahedron 20 30 12 2
Icosahedron 12 30 20 2

นอกจากนี้ polyhedra ยังมีคุณสมบัติ dual กับตัวมันเองหรือ polyhedra แบบอื่นอีกด้วย

Ref:
- http://en.wikipedia.org/wiki/Euler_characteristic
- http://www.light-weaver.com/vortex/5geometry.html


Tags: , , ,

Related posts:

Navier-Stokes Equations ปัญหาที่ยังรอคำตอบ January 31, 2011

Posted by viboon in : Mathematics, Science and engineering , add a comment

1 ใน 7 ปัญหาคณิตศาสตร์ที่ยากและสำคัญที่สุด (The Millennium Prize Problems) ของศตวรรษที่ 21 และยังหาคำตอบในเชิงความสัมพันธ์ไม่ได้คือ Navier-Stokes Equations

คำตอบของสมการนี้จะทำให้เราเข้าใจและสามารถทำนายการปั่นป่วนของของไหลได้อย่างแม่นยำขึ้น สมการนี้ถูกคิดขึ้นในศตวรรษที่ 19 หากแต่ความเข้าใจต่อตัวสมการยังมีอยู่น้อยมาก เพื่อเป็นการพัฒนาวงการวิทยาศาสตร์และคณิตศาสตร์ ผู้ที่สามารถไขความลับของสมการนี้ได้ ทาง Clay Mathematics Institute เสนอเงินรางวัล $1,000,000 และคุณจะมีชื่ออยู่ในตำรา Fluid Mechanics ไปตลอดกาล

ดูปัญหาอื่นๆ ที่ยังคงรอคำตอบได้ที่
http://www.claymath.org/millennium/


Tags: , , , , , ,

Related posts:

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:

Numerical methods: Simpson's Rule October 8, 2010

Posted by viboon in : Mathematics , add a comment

For approximating the integral of f(x)

y = int(f(x), a, b)

%n is even number of segments
d = (b-a)/n;
z1 = 0;
x = a + d;
for i = 1:2:n-1
   z1 = z1 + f(x);
   x = x + 2*d;
end
z2 = 0;
x = a + 2*d;
for i = 2:2:n-2
   z2 = z2 + f(x);
   x = x + 2*d;
end
y = (f(a) + f(b) + 4*z1 + 2*z2)*d/3;

Tags: , ,

Related posts:

Numerical methods: Root finding October 8, 2010

Posted by viboon in : Mathematics , add a comment

Bisection method

if f(a) <= 0
   lo = a; hi = b;
else
   lo = b; hi = a;
end
mid = lo + (hi-lo)/2
while ((mid ~= lo) && (mid ~= hi))
   if f(mid) <= 0
      lo = mid;
   else
      hi = mid;
   end
   mid = lo + (hi-lo)/2
end

Newton's method (a.k.a. Newton-Raphson method)

while err > errthreshold
   dx = f(x)/df(x);
   x = x - dx;
   err = abs(dx/x)*100;
end

Secant method (a finite difference approximation of Newton's method)

while err > errthreshold
   dx = (f(x1)-f(x0))/(df(x1)-df(x0));
   x0 = x1;
   x1 = x1 - df(x1)*dx;
   err = abs(df(x1)*dx/x1)*100;
end

Tags: , , , , , ,

Related posts:

การหา inverse matrix ด้วย LU Factorization และ MATLAB September 25, 2010

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

กำหนดให้ [A][x]=[B] โดยที่ [A] เป็น matrix ขนาด nxn และ [B] ขนาด nx1 โดยตัวอย่างนี้ เราต้องการหา matrix x ที่มีขนาด nx1

ปกติแล้วเรามักคุ้นเคยกับการหา A-1 ซึ่งก็หาได้จาก

A-1 = (1/det(A))*adj(A)T

หรือใน MATLAB จะใช้คำสั่ง

1
inv(A)

หรือ

1
A^(-1)

ก็ได้

แต่การหา x โดยการใช้วิธีนี้ มักมี error เกิดขึ้นอันเนื่องมาจากการปัดเศษ โดยเฉพาะในการคำนวณในลักษณะ numerical computation ซึ่ง error ที่เกิดขึ้นจะสะสม ทำให้คำตอบสุดท้ายผิดเพี้ยนไปได้ ซึ่งหลักการนี้สามารถเห็นได้จาก numerical methods ประเภท implicit ทั่วไป เช่น backward euler หรือ Crank-Nicolson เป็นต้น

วิธีที่จะนำเสนอคือ LU Factorization ซึ่งถือว่าเป็น numerical method ที่เก่าแก่วิธีการหนึ่ง เกิดจากการคิดค้นวิจัยในช่วงปี 1955 ถึง 1965 โดยชื่อวิธีการตั้งตาม C.F. Gauss

การทำ LU Factorization ซึ่งใช้หลักการ Gauss elimination จะแยกการคำนวณออกเป็น 2 ส่วนคือ Forward substitution และ Backward substitution ตามลำดับ

เริ่มจาก กำหนดให้ [A]=[L][U] โดยที่ [L] คือ lower matrix และ [U] คือ upper matrix ซึ่งเมื่อนำเข้าไปแทนในสมการตั้งต้นจะได้ [L][U][x]=[B]

สิ่งที่เราต้องการหาก็คือ [L] และ [U] ซึ่งหาได้จากการใช้ row operation จัดรูป matrix A ให้อยู่ในรูปของ upper matrix ก่อน
[A|I] -> [U|R]
โดยที่ [U]=[R][A]
จากสมการนี้ สังเกตได้ว่า [R] ก็คือ [L]-1 นั่นเอง หรือ [U]=[L]-1[A]

ในขั้นตอนถัดไปคือการทำ Forward substitution โดยกำหนดให้ [z]=[U][x] ดังนั้น เราจะได้ [L][z]=[B] ซึ่งเราสามารถทำการแก้สมการเพื่อหา [z] ได้โดย [z]=[L]-1[B] หรือก็คือ [z]=[R][B]

ต่อมาเป็นการทำ Backward substitution โดยแทน [z] เข้าไปในสมการก่อนหน้า เพื่อใช้ในการหา [x] โดย [U][x]=[z] ขั้นตอนนี้จะต้องหา inverse U ซึ่งเราสามารถหาได้ง่ายเพราะเป็น triangular matrix โดยที่ [x]=[U]-1[z]

นอกจากจะใช้ row operation ดังที่แสดงไว้ก่อนหน้านี้ การหา [U] และ [L] เรายังสามารถใช้การหาค่าใน matrix ได้โดยตรงหากเป็น matrix ขนาดไม่ใหญ่มาก เช่น

สำหรับการใช้ MATLAB เราสามารถหา [L] และ [U] ได้จาก

1
[L,U]=lu(A)

และหากต้องการหา [x] ก็ใช้คำสั่ง

1
x=inv(U)*inv(L)*B

นอกจากนี้ MATLAB ยังมี operator ที่คิดขึ้นมาเฉพาะ MATLAB เพื่อจัดกการกับ error ที่เกิดขึ้นในการหา inverse แบบ inv(A) คือ การใช้ backward slash “\”

หากพิจารณาปัญหาตามตัวอย่างที่กล่าวมาแล้ว [x] สามารถหาได้ด้วยคำสั่ง

1
x=A\B

error ที่เกิดจากการปัดเศษ หรือ roundoff error ที่เกิดขึ้นในขั้นตอนการหา inverse matrix แบบดั้งเดิมนั้นอาจดูไม่เป็นเรื่องสำคัญ หากเป็นการแก้ปัญหาที่ไม่ซับซ้อนมาก แต่สำหรับการคำนวณที่ต้องการความถูกต้องมากๆ อย่างเช่นการแก้ปัญหาด้วย numerical methods แบบต่างๆ เช่น finite difference method, finite element method หรือ finite volume method การสะสมของ roundoff error นำไปสู่การประมาณค่าสุดท้ายที่ผิดพลาดอย่างมากได้

รูปตัวอย่างจาก: http://www.alkires.com/103/6_lu.pdf


Tags: , , , , , ,

Related posts: