Monday, November 18, 2013

Belajar optimisasi di MATLAB

Sebenarnya untuk kreatif g perlu mahal, asalkan kita banyak baca. Berikut ini saya akan membahas sebuah contoh dasar pada buku kalkulus James stewart edisi 5 pada bab 4.7 yakni soal latihan nomor 57.


Yang hasilnya saya gambar ulang menjadi


Pada gambar di atas, terlihat jelas bahwa sudut \(\theta\) dapat dinyatakan sebagai fungsi \(x\) yakni $$ \begin{equation} \theta = 180^o -tan^{-1} \left ( \frac{2}{x} \right ) - \tan^{-1} \left( \frac{5}{3-x} \right ) \end{equation} $$ yang nilai maksimumnya dapat dinyatakan sebagai $$ \begin{eqnarray} \frac{d \theta}{dx } &= & 0 + \frac{2}{x^2} \frac{1}{1 + \frac{4}{x^2} } - \frac{5}{(3 - x)^2} \cdot \frac{1}{1 + \frac{25}{(3 - x)^2}} \\ & = & \frac{2}{x^2 + 4} - \frac{5}{(3 - x)^2 + 25} = 0 \end{eqnarray} $$ atau \begin{eqnarray} 3 x^2 + 12 x - 12 = 0 \end{eqnarray} yang akarnya adalah \(x = - 6 \pm 6\sqrt{2} \). Di mana yang memenuhi adalah \(x = -6 + 6 \sqrt{2} = 6 (\sqrt{2} - 1) \approx 2.4\) Untuk mengilustrasikan jawaban tersebut, saya sudah membuat sebuah program sederhana dalam bahasa MATLAB yang sourcenya pembaca bisa coba di rumah, ditulis dengan menggunakan MATLAB 2009. Adapun sourcenya adalah:
% jawaban soal nomor 57 BAB 4.7 dari buku James Stewart kalkulus
% edisi 5
% author: Mohammad Fajar
function soal57
clc; 
f = figure('menubar', 'none', 'resize', 'off', ... 
    'position', [200, 100, 800, 450]); 
ax1 = axes('units', 'pix', 'position', [50 50 340 270]... 
   ,  'xlim', [0 2*pi], 'ylim', [-1.5, 1.5]); 

ax2 = axes('units', 'pix', 'position', [430 50 340 270]... 
   ,  'xlim', [0 2*pi], 'ylim', [-1.5, 1.5]); 

slider =  uicontrol('style', 'slider', 'position', [50 380 90 20]); 
text_ = uicontrol('style', 'text', 'position', [150 380 40 20 ]); 

set(slider, 'value', 0.1);
set(slider, 'sliderstep', [0.05 0.05]);
ss = 0:0.01:3; 
theta = 180 - atand(2./ss) - atand(5./(3- ss));
maksimum = max(theta); % jawabannya
while ishandle(f)
    xx = get(slider , 'value');
    xx = xx * 3;
    set(text_, 'string', num2str(xx)); 
    trigx = [0 0;5 2;0 0]; 
    trigy = [0 3-xx;0 3;3-xx 3];
    z = [1 1;1 1;1 1];
    p = patch(trigx, trigy, z, 'parent', ax2); 
    axis([-1 ,6, -1, 4]);
    thetax = 180 - atand(2./xx) - atand(5./(3 - xx)); 
    p2 = plot(ss, theta, 'parent', ax1);
    xlabel(ax1, 'x');
    ylabel(ax1, '\theta'); 
    hold on; 
    p3 = line([xx xx], [0 thetax], 'linewidth', 4, 'parent', ax1);
    tt = text(xx+.1,thetax-3,['\theta = ', num2str(thetax)],'parent', ax1 ...
        , 'fontsize', 9);
    tt2 = text(0.8,3-xx,['\leftarrow \theta = ', num2str(thetax), '^o'],'parent', ax2 ...
        , 'fontsize', 12);
    % ini hanya untuk mengakali karena step slider dan step linspace pada 
    % MATLAB tidak sama :)
    if floor(thetax) == floor(maksimum)
        set(p3, 'color', 'r');
        set(p, 'facecolor', 'y'); 
    end
    drawnow;
    if ishandle(p), delete(p);end
    if ishandle(p2), delete(p2);end
    if ishandle(p3), delete(p3);end
    if ishandle(tt), delete(tt);end
    if ishandle(tt2), delete(tt2);end
end
adapaun hasil runningnya adalah:


Wednesday, November 13, 2013

Memahami garis singgung dengan MATLAB

Berikut ini saya berikan contoh program untuk memahami makna filosofis dari garis singgung.
function garis_singgung
% author: mohammad fajar :)
clc; 
f = figure('menubar', 'none', 'resize', 'off'); 
axes('units', 'pix', 'position', [50 50 450 300]... 
   ,  'xlim', [0 2*pi], 'ylim', [-1.5, 1.5]); 
slider =  uicontrol('style', 'slider', 'position', [50 380 90 20]); 
text = uicontrol('style', 'text', 'position', [150 380 40 20 ]); 
t = 0:.1:2*pi;  
y1 = sin(t);  
y2 = cos(t); 
plot(t, y1 , t, y2);
hold on ;
line([0 2*pi], [0 0]);
hold on ;
while ishandle(f)
    value = get(slider , 'value');
    value = value * 2 * pi; 
    set(text, 'string', num2str(value)); 
    hasil = cos(value); 
    yy = (t - value ).* hasil + sin(value) ;
    p = plot(t, yy, 'r' , 'linewidth', 2); 
    p1 = line([value value ] , [0 cos(value)], 'linewidth', 4);
    p2 = plot(value, sin(value),'r.' , 'markersize', 23); 
    % method drawnow tidak mengenali settingan di atas
    axis([0 , 2*pi, -1.5, 1.5 ]); 
    drawnow; 
    if ishandle(p), delete(p); end
    if ishandle(p1), delete(p1); end
    if ishandle(p2), delete(p2); end
end
end
Adapun hasilnya adalah pada gambar berikut:

Friday, November 1, 2013

Belajar monte carlo dengan MATLAB untuk menghitung luas kurva

Berikut ini contoh program yang saya buat pake MATLAB untuk menghitung luas kurva "apa saja" dengan menggunakan metode Monte Carlo.


Yang source codenya 
function test_monte_carlo
clear all ;
clc;
f = figure('menubar', 'none', 'resize', 'off' ); 

ax = axes('units', 'pix', 'position', [40 50 340 340]);

uicontrol('units', 'pix', 'position', [400 360 60 25],... 
    'style', 'text','string','f(x) :' ,  'fontsize', 12, 'horizontalAlignment', 'left', ...
    'fontsize', 12, 'fontweight', 'bold', 'backgroundcolor', get(f, 'color'));

% input fungsi
inputFungsi = uicontrol('units', 'pix', 'position', [450 360 100 25],... 
    'style', 'edit', 'fontsize', 12);

uicontrol('units', 'pix', 'position', [400 330 60 25],... 
    'style', 'text','string','N :' ,  'fontsize', 12, 'horizontalAlignment', 'left', ...
    'fontsize', 12, 'fontweight', 'bold', 'backgroundcolor', get(f, 'color'));

inputN = uicontrol('units', 'pix', 'position', [450 330 100 25],... 
    'style', 'edit', 'fontsize', 12, 'fontweight', 'bold');

% keterangan hasil
ket_hasil  = uicontrol('units', 'pix', 'position', [400 280 100 30],... 
    'style', 'text', 'fontsize', 12, 'fontweight', 'bold', ...
    'string', '0' , 'backgroundcolor', get(f, 'color'));

% hitung button
uicontrol('units', 'pix', 'position', [400 250 100 30],... 
    'style', 'pushbutton', 'fontsize', 12, 'fontweight', 'bold', ...
    'string', 'HITUNG', 'Callback', @fungsi_hitung);

% reset button
uicontrol('units', 'pix', 'position', [400 210 100 30],... 
    'style', 'pushbutton', 'fontsize', 12, 'fontweight', 'bold', ...
    'string', 'RESET', 'Callback', @reset_fungsi);

    function fungsi_hitung(varargin)
        N = get(inputN , 'string'); 
        N = str2double(N); 
        x = 0:.01:2*pi; 
        yy = get(inputFungsi, 'string'); 
        y = eval(yy, x); 
        plot(x,y, 'parent', ax, 'linewidth', 3);
        hold on;
        miny = min(y); 
        maxy = max(y);
        deltax = linspace( 0 , 2*pi ,10000); 
        deltay = linspace( miny,maxy, 10000); 
        a = 0;
        for i=1:N,
            m = randi(length(deltax),1);
            n = randi(length(deltay),1);
            x = deltax(m); 
            y = deltay(n);
            if  y  >= 0 
                if y  <= eval(yy, x)
                    a = a + 1;
                    plot(x,y, 'g*', 'parent', ax);
                    hold on ;
                else
                   plot(x,y, 'r*', 'parent', ax);
                   hold on ;
                end
            else
                if y  >= eval(yy, x)
                    a = a + 1;
                     plot(x,y, 'g*', 'parent', ax);
                     hold on ;
                else
                     plot(x,y, 'r*', 'parent', ax);
                     hold on ;
                end
            end
        end
        axis tight; 
        luas =  (a/N) * ( (maxy- miny) * (2*pi - 0)); 
        set(ket_hasil, 'string', num2str(luas)); 
    end
    
    function reset_fungsi(varargin)
        cla reset; 
        set(ket_hasil, 'string', '0'); 
    end
end