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:

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: