エカつきのブログ = Eka tsuki no blog

Program Matlab untuk persamaan difusi 1 dimensi

clc;
clear;
disp('Program Persamaan Difusi Satu Dimensi Non-Dimensionalisasi')
disp('-------------------------------------')

%% Rekam PDP dengan syarat awal dan syarat batas
flag=0;
while flag==0
    L=input('Masukkan panjang kawat = ');
    if L<=0
        disp('Nilai panjang kawat harus positif')%syarat positif
    else
        break
    end
end

while flag==0
    alpha=input('Masukkan konstanta alpha = ');
    if alpha<=0
        disp('Nilai alpha harus positif')%syarat positif
    else
        break
    end
end

while flag==0
    h=input('Masukkan banyaknya space step = ');
    if h<=0 || ceil(h)-h>0
        disp('Nilai space step (h) harus bilangan bulat positif')%syarat bulat positif
    else
        a=L/h;
        disp(['Panjang segmen dx= ',num2str(a)])
        break
    end
end

while flag==0
    i=input('Masukkan banyaknya time step = ');
    if i<=0 || ceil(i)-i>0
        disp('Nilai time step (i) harus bilangan bulat positif')%syarat bulat positif
    else
        break
    end
end

disp('---Input Syarat Awal---')
fx=input('Masukkan nilai syarat awal= ');

disp('---Input Syarat Batas---')
T1=input('Masukkan suhu batas pertama = ');
T2=input('Masukkan suhu batas kedua = ');

%% Transformasi ke PDP non-dimensional
%non-dimensionalisasi panjang
dx=1/h;
%non-dimensionalisasi waktu 
dt=1/i;
%non-dimensionalisasi suhu
Tm1=T1/max(T1,T2);
Tm2=T2/max(T1,T2);

% Syarat Stabil
S=(alpha*dt)/dx^2;
if S>1/2
    disp('____________________')
    disp('Kondisi tidak stabil')
else
    disp('____________________')
    disp('Kondisi stabil')
end

%% Grid komputasi
%syarat awal
for j=1:h+1
    x(j)=(j-1)*dx;
    u(j,1)=fx;
end
%syarat batas
for n=1:i+1
    u(1,n)=Tm1;
    u(h+1,n)=Tm2;
    waktu(n)=(n-1)*dt;
end

%% Menyelesaikan dengan metode FTCS
for n=1:i
    for j=2:h
        u(j,n+1)=S*u(j+1,n)+S*u(j-1,n)+(1-2*S)*u(j,n);
    end
end
%plot(x, u')
mesh(x,waktu,u')
xlabel('Panjang')
ylabel('Waktu')
zlabel('Suhu')

%% Mengembalikan PDP dimensional
panjang=x'*L;
suhu=u(:,n+1)*max(T1,T2);
disp('-----------------------------')
disp('    Panjang       Suhu')
disp('-----------------------------')
[panjang, suhu]