option nolet print "------------------------------------------" print " One-Dimensional Steady Heat conduction " print " solve: ode: k*d/dx(dT/dx) =0(version 1.0)" print "------------------------------------------" print " rectangular slab ,divide=9 ,no=8 " print " i=1, T=100°C,i=10 ,T=200°C ,k=120 W/cm-K " print "------------------------------------------" print "algo: TDMA = TriDiagonalMatrixAlgorithm " print "peter.vlasschaert@gmail.com (27/06/2015) " print " model (vesion 1.0),sub version " print "------------------------------------------" print print "------------------------------------------" print " a(i)*x(i-1)+b(i)*x(i)+c(i)*x(i+1) = d(i) " print " a(1)=0, c(no)=0 => " print " x(i)=gamma(i)-c(i)*x(i+1)/beta(i) " print "------------------------------------------" print " d(i)=aa(i,j)*x(j) " print "------------------------------------------" print input prompt " Number of equations = ":no input prompt "sel =1, Heat One-Dimensional (conduction) = ":sel if sel=1 then md_b=2 ad_c=-1 bd_a=-1 else input prompt " Main Diagonal = ":md_b ! -2 input prompt " Above Diagonal = ":ad_c ! 1 input prompt " Below Diagonal = ":bd_a ! 1 end if input prompt " Temperature left side slab = ":t1 input prompt " Temperature right side slab = ":t2 dim a(1),b(1),c(1),d(1) mat redim a(no),b(no),c(no),d(no) call in_tridiag(no,a(),b(),c(),md_b,ad_c,bd_a) !=> aa*x=d call in_d(no,d(),t1,t2) !--------------------------- ! solve TDMA !--------------------------- ! x(i) -> i=1..no !--------------------------- dim beta(1),gamma(1) mat redim beta(no),gamma(no) call beta_gamma(no,a(),b(),c(),d(),beta(),gamma()) !///////////////solve (back substitution) ////////////////// dim x(1) mat redim x(no) call back_solve (no,beta(),gamma(),c(),x()) dim aa(1,1) mat redim aa(no,no) print " d(i)=aa(i,j)*x(j) => show: d(i) " print mat print d print " matrix aa(i,j) => " print call show_tridiag(no,a(),b(),c(),aa(,)) print " column temperature(j) => " print mat print x end sub beta_gamma(no,a(),b(),c(),d(),beta(),gamma()) beta(1) = b(1) gamma(1)= d(1)/beta(1) for j=2 to no beta (j)=b(j)-a(j)*c(j-1)/beta(j-1) gamma(j)=(d(j)-a(j)*gamma(j-1))/beta(j) next j end sub sub back_solve (no,beta(),gamma(),c(),x()) x(no)=gamma(no) ! back substitution : solve n->n-1->....1 for k=no-1 to 1 step -1 x(k)=gamma(k)-c(k)*x(k+1)/beta(k) next k end sub sub show_tridiag(no,a(),b(),c(),aa(,)) for i=1 to no !=> main diagonal aa(i,i)=b(i) !=>above diagonal if ibelow diagonal if i>=2 then aa(i,i-1)=a(i) end if next i mat print aa end sub sub in_tridiag(no,a(),b(),c(),md_b,ad_c,bd_a) !=> coef main diag for i=1 to no b(i)=md_b next i !=> coef above diag for i=1 to no-1 c(i)=ad_c next i !=> coef below diag for i=2 to no a(i)=bd_a next i end sub sub in_d(no,d(),t1,t2) d(1)=t1 for i=2 to no-1 d(i)=0 next i d(no)=t2 end sub