update : 26/07/2016
Thiele-Geddes sumproducts
download code: Thielegeddessumproducts
output view : view Thielegeddessumproducts
option nolet
print"**************************************************************"
print" example | ( n=3 ) input: A0=1;A1=2;A2=3;A3=4 "
print"--------- "
print" called : Thiele-Geddes sumproducts "
print"**************************************************************"
print" sum = A0*A1*A2*A3+A1*A2*A3+A2*A3+A3+1=65 "
print"**************************************************************"
print
print
dim a(0 to 1)
dim b(0 to 1)
n=3 ! to change n you can make any computation (with this values above)
mat redim a(0 to n)
mat redim b(0 to n)
!**************************
! input values A0,A1,A2,A3*
!**************************
for i=0 to n
a(i)=i+1
next i
!-----------------------
! products:A0*A1*....*Ai
!-----------------------
prod=1 !init
for i=0 to n
prod=1
for j=i to n
prod=prod*a(j)
next j
b(i)=prod
next i
print "******************"
print "matrix of products"
print "******************"
mat print b
print
print"!---------------------"
print"! (sum of products) +1"
print"!---------------------"
som=1 !init
for i=0 to n
som=som+b(i)
next i
print "sum = ";som
end
graphics : sub routines (below) , general use when draw graphics on the screen.
sub draw_asp(xwin,ywin,asp)
set window -xwin,xwin,-ywin,ywin
! control aspect ratio for circle not ellips
ask pixels px,py
if px>py then
asp=px/py
else
asp=py/px
end if
end sub
SUB draw_axes(xmin,xmax,ymin,ymax,nt)
LET ntick = nt ! number of tick marks
! distance between tick marks on x-axis
LET dx = (xmax - xmin)/ntick
! distance between tick marks on y-axis
LET dy = (ymax - ymin)/ntick
! include margin
SET WINDOW xmin - dx,xmax + dx,ymin - dy,ymax + dy
LET x0 = max(0,xmin)
LET y0 = max(0,ymin)
IF ymin*ymax < 0 then
LET y0 = 0
ELSE
LET y0 = ymin
END IF
PLOT LINES: xmin,y0;xmax,y0 ! horizontal axis
PLOT LINES: x0,ymin;x0,ymax ! vertical axis
LET tx = 0.1*dy ! size of tick mark
LET ty = 0.1*dx
FOR itick = 0 to ntick
LET x = xmin + itick*dx
PLOT LINES: x,y0 - tx; x,y0 + tx ! draw ticks on x axis
LET y = ymin + itick*dy
PLOT LINES: x0 - ty,y; x0 + ty,y ! draw ticks on y axis
NEXT itick
END SUB
sub smallcir_dr(s1,xf,yf,scale,asp,z3$)
set color z3$
for theta=0 to 2*3.14 step 0.01
plot points:s1*cos(theta)+xf,s1*sin(theta)*asp+yf
next theta
end sub
================================================================
The use of the 'mat' command.
1e) sub
sub somprodmat(i,j,m(,),n(,),c(,),bg,be)
som=0
for k=bg to be
som = som + m(i,k)*n(k,j)
next k
c(i,j) = som
end sub
sub somprodvec(i,m(,),z(),cc(),bg,be)
som=0
for k=bg to be
som = som + m(i,k)*z(k)
next k
cc(i) = som
end sub
2e) 'see output view : below ,sub below'
option nolet
print "*****************************"
print " info ,'mat' command "
print " [email protected] "
print "*****************************"
print
dim a(2,1)
a(1,1)= 7
a(2,1)=8
dim b(1,2)
b(1,1)=5
b(1,2)=6
dim c(2,2)
c(1,1)=1
c(1,2)=2
c(2,1)=3
c(2,2)=4
dim e(2,2)
e(1,1)=1
e(1,2)=2
e(2,1)=3
e(2,2)=4
dim d(1,2)
print " given = matrix b,c "
print " b= "
mat print b
print " c= "
mat print c
!*************************************
! multiplication two matrices
! b(1,2)*c(2,2) => d(1,2)
! (1,2)*(2,2) contract to (1,2)
!*************************************
! 1e) (1,2)*(2,2) = (1,2)
!*************************************
call prod_mat(1,2,2,b(,),c(,),d(,))
print " matrix = d =>multiplication matrix 'b*c' "
print " dimension => dim d(1,2) "
mat print d
dim cmat(1,1)
!*************************************
! multiplicatication two matrices
! d(1,2)*a(2,1) = > number
!*************************************
! 2e) (1,2)*(2,1) = (1,1) = number
!*************************************
print " given = matrix a "
print " a= "
mat print a
print "matrix = cmat =>multiplication matrix 'd*a'"
call prod_mat(1,2,1,d(,),a(,),cmat(,))
print "result => matrix multiplication = 'b*c*a' "
print " cmat ' number ' = "
print
mat print cmat
print
kk=1
ll=-4
mat redim cmat(2,2)
i=2
j=2
print " sum of two matrix kk= ";kk;" ll= ";ll
print " given = matrix c,e "
print "c= "
mat print c
print "e= "
mat print e
call som_mat(kk,ll,i,j,c(,),e(,),cmat(,))
print " result cmat = ";kk;"*c+ ";ll;"*e"
print
mat print cmat
end
sub prod_mat(i,k,j,amat(,),bmat(,),cmat(,))
mat redim amat(i,k),bmat(k,j),cmat(i,j)
mat cmat = amat*bmat
end sub
sub som_mat(kk,ll,i,j,amat(,),bmat(,),cmat(,))
!use :kk,ll float
! i,j integer
mat redim amat(i,j),bmat(i,j),cmat(i,j)
mat amat =kk*amat
mat bmat =ll*bmat
mat cmat = amat+bmat
end sub
==================================================================
download code: basicmatrixconstruct
output view :viewbasicmatrixconstruct
special functions ( example : series : Bessel function of order n (integer) = Jn(t))
n = integer include 0,1,2,3,....
t = ?
coef = t^n/(2^n*n!)
=> Jn(t) = coef*(1-t^2/(2*(2*n+2))+t^4/(2*4*(2*n+2)*(2*n+4))-.........)
option nolet
! special functions(Legrendre,Bessel,gamma......"theoretical chemistry=schrodinger eq")
print "******************************************************"
print " chemical engineering "
print " Bessel n th order (integer) version 1.0 "
print " [email protected] , 24-11-2015 "
print "******************************************************"
print
!-------------------------------------------------------------
! ref : schaum's outline series
! laplace transforms !(time => freq domain)
! !( inverse laplace transform ' freq => time domain ')
! Murray R.SPIEGEL
! McGraw-Hill BOOK COMPANY
! including 450 solved problems ( page 7,edition 1965)
! general : every edition can be used
!------------------------------------------------------------
!************************************************************
! 1e)zero order besselfunction. = b0 = f_Bessel0
! 2e) n (integer order) besselfunction = bn = b_n
!************************************************************
! example =>
! serie infinite number of terms (4 terms here => q=4,t= ?)
!************************************************************
t=1
t2=t*t
t4=t2*t2
t6=t4*t2
b0 = 1-t2/4+t4/(2*2*4*4)-t6/(2*2*4*4*6*6)
print " b0 (bessel zero order) = ";b0
!*****************************************************
q=3 ! first term ( som = 1 ) init => q=4
!*****************************************************
! programming series :bessel n(integer) order =>
!*****************************************************
call f_Bessel0(t,q,som)
print " f_Bessel0 (bessel zero order )= ";som
p=0 ! integer
! without coef
p2 = 2*(2*p+2)
p4 = p2*4*(2*p+4)
p6 = p4*6*(2*p+6)
bn = 1-t2/p2+t4/p4-t6/p6
print " bn (integer) = ";bn
call f_Besseln(p,t,q,som)
print " bn (integer) = ";som
pp = 0
call faculteit(pp,pr)
print " pr = ";pr
call coefbesseln(t,pp,pr,coef_besseln)
b_n = som*coef_besseln
print " with coef besseln (bessel n(integer) order = ";b_n
end
sub f_Bessel0(t,q,som)
product=1
som =1 ! first term is 1 otherwise use 0
for m=1 to q
product = 1
for n=1 to m
product = product*(2*n)*(2*n)
next n
som = som + ((-1)^m*t^(2*m))/product
next m
end sub
sub f_Besseln(p,t,q,som)
! without coef
product=1
som =1 ! first term is 1 otherwise use 0
for m=1 to q
product = 1
for n=1 to m
product = product*(2*n)*(2*p+(2*n))
next n
som = som + ((-1)^m*t^(2*m))/product
next m
end sub
sub faculteit(pp,pr)
if pp = 0 then
pr= 1
else
pr = 1
for i=1 to pp
pr=pr*i
next i
end if
end sub
sub coefbesseln(t,pp,pr,coef_besseln)
coef_besseln = t^pp/(2^pp*pr)
end sub
=======================================================
def :kgv=kleinste gemeen veelvoud(= dutch) = lcd( least common divider)
ggd=grootste gemene deler (= dutch) = gcd( greatest common divider)
lemma : ggd(a,b)*kgv(a,b)=a*b
option nolet
declare function ggd,kgv
a=3
b=26
c=a*b
ant1=ggd(a,b)
print "grootste gemene deler(";a;",";b;") = ";ant1
ant2=kgv(ant1,a,b)
print "kleinste gemene veelvoud(";a;",";b;") = ";ant2
! floor are the same int
f=int(2.501)
print f
end
def kgv(d,a,b)
c=a*b
kgv=c/d
! ggd(a,b)*kgv(a,b)=a*b
end function
def ggd(d,e)
do
tmp = mod(d,e)
d = e
e = tmp
loop while e <> 0
ggd = d
end function
! gcd simple (greatest common divisor)
! find the common factors of two numbers
! number m,number n
! print "list: use matrix command"
! find: gcd(m,n)
option nolet
m = 63
n = 14
print "***********************************"
print " peter.vlasschaert ,26/07/2016 "
print " gcd:(g)reatest (c)ommon (d)ivisor "
print "***********************************"
print
print " m = number ";m
print " n = number ";n
print
dim fm(1),fn(1)
!fm=factors of m: append to list fm
n1=0
for i=1 to m
if mod(m,i)=0 then
n1=n1+1
mat redim fm(n1)
fm(n1)=i
end if
next i
print " factors of fm = "
print
mat print fm
!***********************************
! result fm : 1,3,7,9,21,63
!***********************************
!fn=factors of n: append to list fn
m1=0
for i=1 to n
if mod(n,i)=0 then
m1=m1+1
mat redim fn(m1)
fn(m1)=i
end if
next i
print " factors of fn = "
print
mat print fn
!***********************************
! result fm : 1,2,7,14
!***********************************
dim cf(1)
! cf=common factors of " fm 'and' fn "
! at least one common divisor '1'
l1=0
n_l1=size(fn)
m_l1=size(fm)
zn = min(m_l1,n_l1)
for i=1 to zn
if fm(i)=fn(i) then
l1=l1+1
mat redim cf(l1)
! cf(l1) = fn(i),choose this one is also possible
cf(l1)=fm(i)
end if
next i
print "list 'cf' of common factors of 'fm and fn' "
print
mat print cf
!*****************
! result cf:1,7
!*****************
!
! find highest common factor (last in the list)
!
n_cf = size(cf)
gcd = cf(n_cf)
print " gcd(";m;",";n;")= ";gcd
!********************
! result gcd(63,14)=7
!********************
end
recursion (Wikipedia),factorial,sinus
---------------------------------------------------------
see above: other way to get factorial definition.
' sub faculteit(pp,pr) '
!*****************************************************
! recursion : factorial function
! fac(n) = 1*2*3...(n-1)*n
! recursion : sinus function ? How see below
! sin(x) = x - x^3/fac(3)+x^5/fac(5)-x^7/fac(7)+...
! recursion : cosinus function ? How see below
! cos(x) = 1 - x^2/fac(2)+x^4/fac(4)-x^6/fac(6)+
!*****************************************************
print "*********************************************"
print " recursion : sin function ,cos function 1.1 "
print " [email protected] , 17/05/2016 "
print " ********************************************"
print
option nolet
declare function fac,a,b
print "fac(";0;")=";fac(0)
for k = 0 to 3
ii = 2*k+1
print "fac(";ii;")=";fac(ii)
next k
print
print " sin(pi/4) = ";sin(pi/4)
print " sin(3.14159/4) = "; sin(3.14159/4)
x = 3.14159/4
print
print " x = ";x
print
nu = 4 ! first term = x;second term = - x^3/fac(3),.....
print " number of terms for 'recursion' = ";nu
print "*************************************************"
print " sin(x) = x - x^3/fac(3)+x^5/fac(5)-x^7/fac(7) *"
print "*************************************************"
print
print " sin(pi/4) = ";x-x^3/fac(3)+x^5/fac(5)-x^7/fac(7)
print "********************"
print " recursion : sin(x) "
print "********************"
print
print "*************************************************"
print " t(k) = (-1)^(k-1)*x^(2*k-1)/fac(2*k-1) "
print " next term : k -> k+1 "
print " t(k+1) = (-1)^k*x^(2*k+1)/fac(2*k+1) "
print "*************************************************"
print " t(k+1) = a(k,x) * t(k) "
print " ************************************************"
print " a(k,x) = t(k+1)/t(k) = - x^2/((2*k)*(2*k+1)) "
print "*************************************************"
sum_sin = x ! first term 'sinus serie '
term = x ! first term
for k = 1 to nu-1 ! nu-1=k => t(k+1)=t(nu) terms
! recursion ' inside '
term = a(x,k)*term
sum_sin = sum_sin + term
next k
print
print "recursion : sin(";x;") = ";sum_sin
print " sin(";x;") = sqr(2)/2 = ";sqr(2)/2
print
for k = 1 to 3
ii = 2*k
print "fac(";ii;")=";fac(ii)
next k
print
print "***********************************************"
print " cos(x) = 1 - x^2/fac(2)+x^4/fac(4)-x^6/fac(6) "
print "***********************************************"
print " (sum_cos)^2+(sum_sin)^2 = 1 'identity' "
print "***********************************************"
sum_cos = sqr(1-(sum_sin)^2)
print " cos(";x;") = ";sum_cos
print " serie cos(";x;")= ";1 - x^2/fac(2)+x^4/fac(4)-x^6/fac(6)
sum_cos = 0 ! classic ' sum = 0 '
term = 1 ! first term
for k = 1 to nu-1 ! nu-1=k => t(k+1)=t(nu) terms
! recursion ' inside '
sum_cos = sum_cos + term
term = b(x,k)*term
next k
print "********************"
print " recursion : cos(x) "
print "********************"
print
print " ************************************************"
print " b(k,x) = t(k+1)/t(k) = - x^2/((2*k)*(2*k-1)) "
print "*************************************************"
print " recursion : cos(";x;") = ";sum_cos
print " cos(";x;") = sqr(2)/2 = ";sqr(2)/2
end
def a(x,k)
a = - x^2/((2*k)*(2*k+1))
end def
def b(x,k)
b = - x^2/((2*k)*(2*k-1))
end def
def fac(i)
!***********************************************************
! recursion ' inside '
! https://en.wikipedia.org/wiki/Recursion_(computer_science)
!***********************************************************
if i = 0 then
fac = 1
else
fac = i*fac(i-1)
end if
end def
output view : view recursionsincosfunc1point1
TableForm : mat (usefull for finite elements )
-----------------------------------------------------------------------
! put smallmatrix in bigmatrix (for every element)*
! very handy tool for 2d matrix (table form) *
!**************************************************
! [email protected] , 30/05/2016 *
!**************************************************
print "**********************"
print "* tableform : matrix *"
print "**********************"
print " n = k1*(i-1) + j "
print "*************************"
print " example: k1 =5 "
print "*************************"
print " i,j =counter big matrix "
print "*************************"
!*******************************************************
! why : you put a other matrix in column 8 and then make
! the sum of every element from column 7 and column 8
! and put the sum in column 9 (you can do it )
!*******************************************************
option nolet
pos1 = 3
pos2 = 1
pos3 = 4
print " position => "
print "position(1)= "; pos1
print "position(2)= "; pos2
print "position(3)= "; pos3
print " matrix a => "
dim a(3,3)
!values matrix a ( 'read - data' statement)
for i=1 to 3
for j=1 to 3
read a(i,j)
next j
next i
mat print a
! dim ' big matrix '
k1 = 5
k2 = 5
k3 = k1*k2
dim k(1,1)
mat redim k(k1,k2)
mat k = zero
dim table(1,1)
mat redim table(k3,9)
! column 1
for i=1 to k3
table(i,1) = i
next i
! column 3,4
for i=0 to k1-1
for j=1 to k1
table(j+i*k1,3)=i+1
table(j+i*k1,4)=j
next j
next i
! algo
for i=1 to k3
if i=pos1 or i=pos2 or i=pos3 then !test: row ( position mat :a)
for j=1 to k3
if j=pos1 or j=pos2 or j=pos3 then !test: column ( position mat :a)
table(k1*(i-1)+j,5) = i ! column 5
table(k1*(i-1)+j,6) = j ! column 6
table(k1*(i-1)+j,2) = k1*(i-1)+j ! column 2
!*******************************************************************************
!use: example simplex 2d in a bigger matrix (finite element method:flow problem)
! : triangles ( simplex 2d),three points.
!*******************************************************************************
s1=a(1,1)
!*(pos1,pos1)*
k(pos1,pos1)=s1
table(k1*(pos1-1)+pos1,7) = s1 !column 7
s2=a(1,2)
!*(pos1,pos2)*
k(pos1,pos2)=s2
table(k1*(pos1-1)+pos2,7) = s2 !column 7
s3=a(1,3)
!*(pos1,pos3)*
k(pos1,pos3)=s3
table(k1*(pos1-1)+pos3,7) = s3 !column 7
s4=a(2,1)
!*(pos2,pos1)*
k(pos2,pos1)=s4
table(k1*(pos2-1)+pos1,7) = s4 !column 7
s5=a(2,2)
!*(pos2,pos2)*
k(pos2,pos2)=s5
table(k1*(pos2-1)+pos2,7) = s5 !column 7
s6=a(2,3)
!*(pos2,pos3)*
k(pos1,pos3)=s6
table(k1*(pos2-1)+pos3,7) = s6 !column 7
s7=a(3,1)
!*(pos3,pos1)*
k(pos3,pos1)=s7
table(k1*(pos3-1)+pos1,7) = s7 !column 7
s8=a(3,2)
!*(pos3,pos2)*
k(pos3,pos2)=s8
table(k1*(pos3-1)+pos2,7) = s8 !column 7
s9=a(3,3)
!*(pos3,pos3)*
k(pos3,pos3)=s9
table(k1*(pos3-1)+pos3,7) = s9 !column 7
end if
next j
end if
next i
print " table => "
print "number" ,"number"," i"," j","position i","position j","value:matrix a"
print
print 1,2,3,4,5,6,7,8,9,"column"
print
mat print table
data .3,.14,99
data .5,.8,.6
data .1,.11,.3
end
output view : tableformmat
TableForm : mat (relation : between 'index and number' )
-----------------------------------------------------------------------------------------
! find relation between nu=number and index i,j*
!***********************************************
! [email protected], 30/05/2016
!***********************************************
option nolet
print "******************************************************"
print " tableform : relation between ' index i,j and number' "
print "******************************************************"
maxj = 3 ! max length of the block
dim table(1,1)
k=maxj*maxj
mat redim table(k,4) ! make a table
mat table = zero ! make all entries of mat table zero
for i=1 to k
table(i,1) = i ! column '1' number 'nu'
next i
for i=0 to maxj-1
for j=1 to maxj
table(j+i*maxj,2)=i+1 ! column '2' index 'i'
table(j+i*maxj,3)=j ! column '3' index 'j'
next j
next i
print "display table : lengthblock = maxj = ";maxj
print
print " number = nu ","index 'i'","index 'j'"
mat print table
print
print "*************************************"
print "*************************************"
print "* nu = maxj*(i-1)+j Formula (1)"
print "*************************************"
print "*************************************"
print
ii = 3
print " given = i => ";ii
jj = 2
print " given = j => ";jj
nu1 = maxj*(ii-1)+jj
print " use Formula (1) 'see table to veryfy' "
print " find = nu => ";nu1
print
print "**************************************"
print "**************************************"
print " find inversion for 'Formula (1)' "
print "**************************************"
print " algo : 'Formula (2)' "
print " index 'j' => "
print " j = nu mod maxj "
print " if j = 0 then j = maxj "
print " if j <> 0 then j = nu mod maxj "
print " index 'i' => "
print " i = (nu-j)/maxj + 1 "
print "**************************************"
print "**************************************"
print
nu2 = 9
print " given = nu => ";nu2
!******************************
! implementation of Formula (2)
!******************************
j = mod(nu2,maxj)
print " result j 'check' j=0 or j<>0 => ";j
if j = 0 then
j = maxj
else
j=mod(nu2,maxj)
end if
print
print " check table to veryfy => "
print
print " find = j => ";j
i = (nu2-j)/maxj + 1
print " find = i => ";i
end
output view : tablerelationnumberindex
How to implement BC ( = Dirichlet conditions ) on a Stifness Matrix.
--------------------------------------------------------------------------------------------------------
! how to implement boundary conditions (stifness matrix = 'k')
! => constant values ' start square system '
! first we have system of equations ( Ax=b)
! => finite elements : use A = 'k' (stifness matrix )
! : use b = 'p' (load vector)
! : use x = 'phi' (physical quantity )
! stifness matrix ' read - data '
!************************************************************
option nolet
! filename = boundaryconditiondirichletversion1.tru
print "*********************************************"
print "* how to implement BC (Dirichlet conditions)*"
print "* on a system of equations version 1.0 *"
print "*********************************************"
print "* [email protected],06/06/2016 *"
print "*********************************************"
read nu ! number of equations
dim k(1,1)
mat redim k(nu,nu)
for i=1 to nu
for j=1 to nu
read k(i,j)
next j
next i
print " stifness matrix => read - data "
mat print k
dim p(1,1)
mat redim p(nu,1)
print " load vector => read - data "
for i=1 to nu
read p(i,1)
next i
mat print p
read nc ! number of constraint
dim phi(1)
mat redim phi(nu) ! phi solution vector ( with constraint)
dim tt(1),tt_val(1)
mat redim tt(nc),tt_val(nc)
for i = 1 to nc
read tt(i)
next i
for i = 1 to nc
read tt_val(i)
next i
!************************************
!*constraint can put solution vector*
!************************************
for i=1 to nc
phi(tt(i))=tt_val(i)
next i
print "-------------------------"
print " before solve the system "
print "-------------------------"
mat print phi
! reduce system : ne = nu - nc
ne = nu-nc ! constraint system
dim k_co(1,1)
mat redim k_co(ne,nu)
dim un_tt(1)
mat redim un_tt(nu)
for i=1 to nu
for ii=1 to ne
!**** tt = {1,4,5}********
if tt(ii) = i then
un_tt(i) = tt(ii)
else
exit if
end if
next ii
next i
!******************************
!******************************
! construct ' data-file '
!**** un_tt = {1,0,0,4,5,0}****
!******************************
!******************************
!*************************************
! constraint stifness matrix = k_co
!*************************************
mat print un_tt
ii=0
for i=1 to nu
if un_tt(i) = 0 then ! we need the zero entries
ii = ii+1 ! values must be put in different rows
for j=1 to nu ! ( all values for the columns)
! read dimension (6 by 6) to a (ne by 6)
k_co(ii,j)=k(i,j)
next j
end if
next i
print " constrained stifness matrix => "
mat print k_co
!*************************************
! constraint load vector = p_co
!*************************************
dim p_co(1,1)
mat redim p_co(ne,1)
ii=0
for i=1 to nu
if un_tt(i) = 0 then ! we need the zero entries
ii = ii+1 ! values must be put in different rows
for j=1 to 1 !( j=1 for load-vector)
! read dimension (6 by 6) to a (ne by 6)
p_co(ii,j)=p(i,j)
next j
end if
next i
print " constrained load vector => "
mat print p_co
print
!*************************************************
! system to solve = 'k_sys '
!*************************************************
dim k_sys(1,1)
mat redim k_sys(ne,nc)
s=0
for i=1 to nu
if un_tt(i) = 0 then
s=s+1
for ii=1 to ne
k_sys(ii,s)=k_co(ii,i)
! k_sys(ii,2)=k_co(ii,3)
! k_sys(ii,3)=k_co(ii,6)
next ii
else
exit if
end if
next i
print " only system to solve : after put in the boundary conditions "
print "k_sys => "
mat print k_sys
!*************************************************
! known values,put opposit (-) sign to load vector
!*************************************************
dim p_remove(1,1)
mat redim p_remove(ne,1)
for ii=1 to nc
for i=1 to ne
p_remove(i,1)=k_co(i,tt(ii))
p_remove(i,1)= - phi(tt(ii))*p_remove(i,1)
next i
! mat print p_remove
next ii
! make a zero vector ' p_zero ' contain = sum of vectors
dim p_zero(1,1)
mat redim p_zero(ne,1)
!mat p_zero = 0
!****************************
! we can use also direct p_co
!****************************
mat p_zero = p_co
print
for ii=1 to nc
for i=1 to ne
p_remove(i,1)=k_co(i,tt(ii))
p_remove(i,1)= - phi(tt(ii))*p_remove(i,1)
next i
mat p_zero = p_zero + p_remove
next ii
!**********************************
! result pp of system we need solve
!**********************************
mat p_co = p_zero
print " load vector for 'k_sys' "
mat print p_co
!solve ( example 'data - file ' phi(2),phi(3),phi(6)'
dim k_inv(1,1)
mat redim k_inv(ne,nu-nc)
mat k_inv = inv(k_sys)
! mat print k_inv
dim result(1,1)
mat redim result(ne,1)
mat result = k_inv*p_co
! example :given 'phi(1),phi(4),phi(5) '
print " solve k_sys => "
mat print result
print "**************************************"
print " solution vector ph(1) to ph(";nu;") "
print " 'k'=stifness matrix 'p'=load vector "
print "**************************************"
ii = 0
for i=1 to nu
if un_tt(i) = 0 then
ii = ii+1
phi(i)=result(ii,1)
else
exit if
! print ii,i,phi(i)
end if
next i
mat print phi
data 6 ! nu=number of equations
! data : stifness matrix ( example 6 by 6)
data 1000,10,20,30,40,50
data 10,2000,21,31,41,51
data 20,21,3000,32,42,52
data 30,31,32,4000,43,53
data 40,41,42,43,5000,54
data 50,51,52,53,54,6000
! load vector
data 100,110,120,130,140,150
! constraint
data 3 ! nc = number of constraint
data 1,4,5 ! matrix tt ( i from phi )
data 5,-7,0 ! matrix tt_val ( ph(i) = values )
end
output view : boundaryconditiondirichletversion1
Parse a string : a$ ( with values inside and some delimiter = '|,;')
---------------------------------------------------------------------------------------------------
'SUB version ' = parse string a$
download code: DELIMITERSTRINGFINALSUBFINAL
output view : parsestringsubfinal
option nolet
!***************************************
! [email protected],08/06/2016
!***************************************
! parse a string : a$ !
!***************************************
print
print "*************************"
print "*** parse a string a$ ***"
print "*************************"
print "****************************************"
print "*[email protected],08/06/2016*"
print "****************************************"
print
a$ = "1.54,1.2,3.56"
print " string a$ = "
print
print a$
print
dim tab$(3)
! selection
tab$(1) = "|"
tab$(2) = ","
tab$(3) = ";"
print " find what type delimiter used : ? ' |,; ' "
print
for j=1 to size(tab$)
p1=pos(a$,tab$(j))
if p1 <> 0 then
pa = j
end if
next j
print
print pa ,tab$(pa)
print " find position first delimiter in a$ "
for i=1 to len(a$)
if a$[i:i]=tab$(pa) then
pb = i
exit for ! only find one : the first
end if
next i
print
print pb,a$[pb:pb]
print
print " find all positions for delimiter in a$"
dim pc(1)
ii = 0
! len : used for the length of the string
for i=1 to len(a$)
if a$[i:i]=tab$(pa) then
ii = ii+1
mat redim pc(ii)
pc(ii) = i
! exit for
end if
next i
nu = size(pc,1)
print
print " number of delimiters = ";nu
print
print " positions of all delimiters = "
print
mat print pc
! search for numbers = nu+1,example three strings
print " strings = "
print
print "first string = "; a$[1:pc(1)-1]
print "second string = "; a$[pc(1)+1:pc(2)-1]
print "thirdth string = "; a$[pc(2)+1:len(a$)]
print
print " convert : strings to values = "
print
print "first number = "; val(a$[1:pc(1)-1])
print "second number = "; val( a$[pc(1)+1:pc(2)-1])
print "thirdth number = "; val( a$[pc(nu)+1:len(a$)])
print
print "***************************************************"
print "*make sum of ";nu+1;" numbers from the string a$ =*"
print "***************************************************"
! init sum
s=0
! first value
s=s+val(a$[1:pc(1)-1])
! between values
for i=1 to nu-1
s=s + val( a$[pc(i)+1:pc(i+1)-1])
next i
! last value
s=s+val( a$[pc(nu)+1:len(a$)])
print
print " the sum of three numbers from string = "
print
print " sum = ";s
end
output view : stringtonumbers
( read -write ) file: use principal 'string parse'
---------------------------------------------------------------------------------------------------
option nolet
!**************************
! data file : testexcel.dat
!**************************
! make this file :
!**************************
! 1,2,3,4
! 5,7,8,9
! 10,11,12,13
!**************************
print "*************************************************"
print " put in path : testexcel.dat "
print "example: C:\Prog\TrueBASIC GOLD v6\testexcel.dat "
print "*************************************************"
print
print " use : for finite elements (data-file) => "
filename$ = "C:\Prog\TrueBASIC GOLD v6\testexcel.dat"
! input prompt "name of file for data : ":filename$
dim a$(1),b$(1)
call read_file(filename$,a$())
call write_file(filename$,a$())
print
mat print a$
!************************************************
! use : above principal 'parse a string'
! example : three strings to parse (a$(i),i=1..3)
!************************************************
print "*********************"
print " content of the file "
print "*********************"
print
print " filename = ";filename$
print
for i=1 to size(a$)
print i;" line row = ";a$(i)
next i
!******************************
! output :
! 1 line row =1,2,3,4
! 2 line row =5,6,7,8
! 3 line row =9,10,11,12
!******************************
end
sub read_file(filename$,array$())
open #4 :name filename$,create newold,org text,access outin
i=0
do while more #4
i=i+1
mat redim array$(i)
line input #4:array$(i)
loop
close #4
end sub
sub write_file(filename$,array$())
open #4 :name filename$,create newold,org text,access outin
erase #4
set #4:margin maxnum
total = Ubound(array$)
for i=1 to total
print #4:array$(i)
next i
close #4
end sub
( write ) file: Use :Ax=b 'matrix and vector' to file
---------------------------------------------------------------------------
see :advanced programma's 67)
option nolet
dim a(1,1),b(1)
!*******************************************
! => Ax = b (system of equations)
!*******************************************
! filename :matrixvectorwritedatafinal.tru
print "************************************"
print " write data : vector = b,matrix = A "
print "************************************"
print " Use : export to other software : "
print " Tecplot ,Matlab,Scilab,Maxima ... "
print "************************************"
print " Use : Ax=b ,concept "
print "************************************"
print
print "*************************************"
print " [email protected],11/06/16"
print "*************************************"
print
mat redim a(5,5),b(5)
! read matrix A( read - data)
for i=1 to size(a,1)
for j=1 to size(a,2)
read a(i,j)
next j
next i
print " Show matrix A = "
print
mat print A
print
print "****************************************"
print " write to : "
print "****************************************"
print
print " ***************************************"
print "* C:\Prog\TrueBASIC GOLD v6\mat_a.dat *"
print "****************************************"
print
filename_mat$ = "C:\Prog\TrueBASIC GOLD v6\mat_a.dat"
call write_matrix(filename_mat$,a(,))
! read vector b (read - data)
print " Show vector b = "
print
for i=1 to Ubound(b)
read b(i)
next i
mat print b
print
print "****************************************"
print " write to : "
print "****************************************"
print
print " ***************************************"
print "* C:\Prog\TrueBASIC GOLD v6\vec_b.dat *"
print "****************************************"
print
filename_vec$ = "C:\Prog\TrueBASIC GOLD v6\vec_b.dat"
call write_vector(filename_vec$,b())
! data Matrix A (read-data)
data 5,2,-4,3,2
data 8,-6,5,1,6
data 4,3,-9,12.5,6
data 8,6,-5,2,8
data 2,3,4,-5,2
! data Vector b (read-data)
data -2,4,3,1,7
end
! lib write data to file(test windows 10)
Sub write_vector(filename$,data())
open #1: name filename$,create newold,org text,access outin
erase #1
set #1:Margin maxnum
mat print #1:data
close #1
end sub
Sub write_matrix(filename$,data(,))
! write matrix to file
open #1: name filename$,create newold,org text,access outin
erase #1
set #1:Margin maxnum
Mat data = TRN(data)
Mat print #1:data
close #1
end sub
( delete ) matrix: delete column and row from matrix a
--------------------------------------------------------------------------------------
option nolet
!****************************************
!****************************************
! delete columns and rows (version 1.0)
! [email protected] ,13/06/2016
!****************************************
!****************************************
! use : together normal and first alone
!****************************************
!****************************************
! select column => col = ? (first )
! select row => row = ? (second)
!****************************************
!****************************************
! deletecolumnandrowmatrixafinal.tru
dim a(1,1)
mat redim a(5,4)
mat read a
! dimension (number rows,number columns)
n = size(a,1)
m = size(a,2)
print "************"
print " matrix A = "
print "************"
print
mat print a
print
print
print "********************************"
print "* DELETE COLUMN AND ROW FROM A *"
print "********************************"
print
col = 2
del_c = col ! delete column 2
print "***********************************"
print "*operation:delete ";del_c;".column*"
print "***********************************"
print
dim v(1,1)
mat redim v(n,m-1)
print "**********************************"
print "* after :operation :col del *"
print "**********************************"
print
call delete_column(n,m,col,a(,),v(,))
mat print v
! v :operation delete row
! v -> w matrix
!**********************************
! ref : matrix a
!**********************************
! n = number rows matrix 'v'
! m-1 = number columns matrix 'v'
! result : delete row = row
! n-1 = number rows matrix 'w'
! m-1 = number columns matrix 'w'
!**********************************
row = 2
del_r =row ! delete row 2
print
print "***********************************"
print "* operation:delete";del_r;".row *"
print "***********************************"
print
dim w(1,1)
mat redim w(n-1,m-1)
print "**********************************"
print "* after :operation :row del *"
print "**********************************"
print
call delete_row(n,m,row,v(,),w(,))
mat print w
! data - mat read (statement)
data 1,2,3,4,5
data 6,7,8,9,10
data 11,12,13,14,15
data 16,17,18,19,20
end
!*******************************************
! sub :first operation then second operation
! a->v then v->w => result : a->w
!*******************************************
sub delete_row(n,m,row,v(,),w(,))
!**********************************
! second: operation 'v->w'
!**********************************
! n = number rows matrix 'v'
! m-1 = number columns matrix 'v'
! result : operation (delete row)
! n-1 = number rows matrix 'w'
! m-1 = number columns matrix 'w'
!**********************************
del_r = row
for j=1 to m-1
for i=1 to n-1
if i >= del_r then
w(i,j) = v(i+1,j)
else
w(i,j) = v(i,j)
end if
next i
next j
end sub
sub delete_column(n,m,col,a(,),v(,))
!**********************************
! first: operation ' a -> v '
!**********************************
! n= number rows matrix 'a'
! m= number columns matrix 'a'
! result : operation (delete column)
! n = number rows matrix 'v'
! m-1 = number columns matrix 'v'
!**********************************
del_c = col ! delete column = col
for i=1 to n
for j=1 to m-1
if j >= del_c then
v(i,j) = a(i,j+1)
else
v(i,j) = a(i,j)
end if
next j
next i
end sub
output view : matrixoperationdel col and row
( delete ) matrix: delete more columns and rows from matrix a
-------------------------------------------------------------------------------------------------
option nolet
!*********************************************
!*********************************************
! delete more columns and rows (version 1.0)
! [email protected] ,14/06/2016
!*********************************************
! moretimesDELETECOLUMNANDROWsimple45forfinal.tru
!*********************************************
!
!*********************************************
!*********************************************
dim a(1,1),col(1),row(1)
mat redim a(4,5),col(3),row(3)
!***************************************
! read - data statement
!***************************************
mat read col
mat read row
mat read a
!***************************************
!***************************************
! dimension (number rows,number columns)
n = size(a,1)
m = size(a,2)
print "************************"
print "* matrix A = *"
print "************************"
print
print "*****************************************"
print " operation mat : delete columns and rows "
print "*****************************************"
print
mat print a
dim v(1,1),w(1,1)
!*****************************
! first step : one column,row
! first run =>
!*****************************
print " first run => "
print
mat redim v(n,m-1)
call delete_column(n,m,col(1),a(,),v(,))
print " delete column = ";col(1)
print
mat print v
mat redim w(n-1,m-1)
call delete_row(n,m,row(1),v(,),w(,))
print " delete row = ";row(1)
mat print w
!***********************
! second thridth combine
!************************
for ii = 2 to size(col)
print ii;"(e the : run "
print
mat redim v(n-(ii-1),m-ii)
call delete_column(n-(ii-1),m-(ii-1),col(ii),w(,),v(,))
print " delete column = ";col(ii)
print
mat print v
mat redim w(n-ii,m-ii)
call delete_row(n-(ii-1),m-(ii-1),row(ii),v(,),w(,))
print " delete row = ";row(ii)
print
mat print w
next ii
!******************************
!****** read - data ***********
!******************************
! data - mat read (statement)*
!******************************
!******************************
!******************************
!column : col
data 2,3,2
! row : row
data 2,3,2
! mat a
data 1,2,3,4,5
data 6,7,8,9,10
data 11,12,13,14,15
data 16,17,18,19,20
end
!*******************************************
! sub :first operation then second operation
! a->v then v->w => result : a->w
!*******************************************
sub delete_row(n,m,row,v(,),w(,))
!**********************************
! second: operation 'v->w'
!**********************************
! n = number rows matrix 'v'
! m-1 = number columns matrix 'v'
! result : operation (delete row)
! n-1 = number rows matrix 'w'
! m-1 = number columns matrix 'w'
!**********************************
del_r = row
for j=1 to m-1
for i=1 to n-1
if i >= del_r then
w(i,j) = v(i+1,j)
else
w(i,j) = v(i,j)
end if
next i
next j
end sub
sub delete_column(n,m,col,a(,),v(,))
!**********************************
! first: operation ' a -> v '
!**********************************
! n= number rows matrix 'a'
! m= number columns matrix 'a'
! result : operation (delete column)
! n = number rows matrix 'v'
! m-1 = number columns matrix 'v'
!**********************************
del_c = col ! delete column = col
for i=1 to n
for j=1 to m-1
if j >= del_c then
v(i,j) = a(i,j+1)
else
v(i,j) = a(i,j)
end if
next j
next i
end sub
output view :moredeletecolrowsforversion1
( delete ) vector: delete elements from vector (sub version)
--------------------------------------------------------------------------------------------
option nolet
!***********************************************
!***********************************************
! delete : elements from vector (version 1.0)
! [email protected] , 15/06/2016
!***********************************************
!filename : MORETIMESDELETEvectorfinalforversionfinal
dim a(1),aa(1)
mat redim a(13),aa(13)
mat read a
print
print "********************************************"
print "* initial ' vector a '(operations,delete) = "
print "********************************************"
print
print " size( vector 'a') = ";size(a)
print
mat print a
dim el(1)
nc = 3 ! number constraint
mat redim el(nc)
mat read el
print "******************************"
print " constraint : delete elements "
print "******************************"
print
for kf = 1 to nc
call del_el_vec(kf,el(),a())
print " delete element from vector 'a' = ";el(kf)
print
mat print a
next kf
print "*****************************************"
print " result = delete elements from vector 'a'"
print "*****************************************"
print
print " size(vector 'a') = ";size(a)-nc
mat redim a(size(a)-nc) ! remove last three positions
print
print " new vector 'a' = "
print
mat print a
!******************************
!****** read - data ***********
!******************************
! data - mat read (statement)*
!******************************
!******************************
!******************************
!
!******************************
!* read vector *
!******************************
data 2,0,0,0,4,0,0,0,2,0,0,0,0
! read : delete el
data 4,8,13
end
sub del_el_vec(kf,el(),a())
for i=1 to el(kf)-2
a(i) = a(i)
next i
for i=el(kf)-1 to size(a)-1
a(i) = a(i+1)
next i
end sub
output view : deleteelementsfromvectorouput
( write ) vector,mat: with delimiter = "," ";" (sub version)
--------------------------------------------------------------------------------------------
option nolet
! write vector to file seperated by 'point comma' = ;
print
print
print "**************************************************"
print "**************************************************"
print "*write : matrix and vector to file with delimiter*"
print "**************************************************"
print "* [email protected] , 18/06/2016 *"
print "**************************************************"
print "**************************************************"
print
print
dim aa(1)
mat redim aa(4)
mat read aa
delimiter$ = ";"
a$ = "C:\Prog\TrueBASIC GOLD v6\testaa.dat"
print "*************************"
print "*write vector:( to file)*"
print "*************************"
print
print " write to file (string) = ";a$
print " delimiter (string) = ";delimiter$
print
print " vector write to file (screen) = "
print
mat print aa
call write_vec_delimiter(a$,delimiter$,aa())
print
print " ************************************"
print " testaa.dat : "
print " "
print " 1;2;3;4 "
print "*************************************"
print
print
!
!*****************
! read-data vector
!*****************
! write matrix to file seperated by 'comma' = ,
data 1,2,3,4
dim bb(1,1)
mat redim bb(3,4)
mat read bb
delimiter$ = ","
b$ = "C:\Prog\TrueBASIC GOLD v6\testbb.dat"
print "*************************"
print "*write matrix:( to file)*"
print "*************************"
print
print " write to file (string) = ";b$
print " delimiter (string) = ";delimiter$
print
print " matrix write to file (screen) = "
print
!**************
!very important
!**************
mat bb = TRN(bb) ! transpose
!**************
mat print bb
print
call write_mat_delimiter(b$,delimiter$,bb(,))
print " ************************************"
print " testbb.dat : "
print " "
print " 1,2,3,4 "
print " 5,6,7,8 "
print " 9,10,11,12 "
print "*************************************"
!******************
! read-data matrix
!******************
data 1,2,3,4
data 5,6,7,8
data 9,10,11,12
end
!******************************************
! lib write : vec and matrix with delimiter
!******************************************
! 1e) write : vector with delimiter
sub write_vec_delimiter(a$,delimiter$,a())
open #1: name a$,create newold,org text,access outin
erase #1
set #1:Margin maxnum
for i=1 to Ubound(a)-1
print #1:a(i);delimiter$;
next i
for i = Ubound(a) to Ubound(a)
print #1:a(i);
next i
close #1
end sub
! 2e) write : matrix with delimiter
sub write_mat_delimiter(b$,delimiter$,bb(,))
open #2: name b$,create newold,org text,access outin
erase #2
set #2:Margin maxnum
n1 = size(bb,1)
n2 = size(bb,2)
!***********************************
! reverse direction of reading
! => We use transpose bb
! normal : first = i other for j
!***********************************
! last element : i control by i loop
!***********************************
for j=1 to n2
for i=1 to n1
if i<>n1 then
print #2:bb(i,j);delimiter$;
else
print #2:bb(i,j);
end if
next i
!*************************************************
! statement :
! very important otherwise :every element one line
!*************************************************
print #2:
next j
close #2
end sub
datafile vector :TESTAA.dat
datafile matrix :TESTBB.dat
output view : writevecandmatwithdelimiter
( delete ) matrix: delete same columns and rows from matrix a
=> split screen.
-------------------------------------------------------------------------------------------------
option nolet
! delete columns from matrix
dim a(1,1),aa(1,1)
mat redim a(5,5),aa(5,5)
!*****************************
! split screen two equal parts
!*****************************
! filename :
!DELETECOLUMNSMATMATRIX55finalfinalrowfinalsplitscreen.tru
open #1: screen 0,0.5,0,1
open #2: screen .5,1,0,1
!*****************************
! open left part of the screen
!*****************************
window #1
print "****************************************"
print " delete :columns and rows => matrix 'a' "
print " [email protected],25/06/2016 "
print "****************************************"
set zonewidth 20
!********************
! read data
!********************
mat read a
!*******************
! initial matrix a
!*******************
print " matrix a = "
mat print a
!************************
!************************
! column deleting matrix
!************************
dim col_n(1)
mat redim col_n(3)
col_n(1) = 2
col_n(2) = 4
col_n(3) = 5
col = size(col_n)
nx = size(a,1)
ny = size(a,2)
print "nx = ";nx
print "ny = ";ny
print
print " deleting positions of columns "
print
mat print col_n
!***********************
!***********************
!***********************
mat aa = zero
mat print aa
dim zz(1,1)
!****************************************************************
! 'first column position'
! delete column = 2 => col_n(1)-0
!****************************************************************
call del_col_matrix(a(,),aa(,),col_n(1)+0,1)
!mat print aa
call resize_mat(zz(,),aa(,),nx,ny)
print " result :: delete column 'a' = ";col_n(1)
print
mat print zz
!****************************************************************
!combine = 'second column position' and 'thridth column position'
!****************************************************************
for i1=2 to col
call del_col_matrix(zz(,),aa(,),col_n(i1)-(i1-1),1)
mat print aa
call resize_mat(zz(,),aa(,),nx,ny-i1)
print " result :: delete column 'a' = ";col_n(i1)
print
mat print zz
next i1
!mat print zz
dim row(1,1)
mat redim row(size(zz,2),size(zz,1))
! reason :same operations on the rows
mat row = trn(zz)
!******************************
! open right part of the screen
!******************************
window #2
set zonewidth 20
! mat 'a' equivalence with mat 'row'
print " matrix 'a' for rows (columm definition) "
print
print " mat row = transpose(zz) "
print
mat print row
dim aaa(1,1),zzz(1,1)
mat redim aaa(size(row,1),size(row,2))
nx=size(row,1)
ny=size(row,2)
print "nx = ";nx
print "ny = ";ny
print
print " deleting :same row positions "
print
dim row_n(1)
mat redim row_n(col)
mat row_n = col_n
mat print row_n
! a:part one
call del_col_matrix(row(,),aaa(,),row_n(1),1)
print
mat print aaa
! b:part one
call resize_mat(zzz(,),aaa(,),nx,ny)
print " result :: delete column 'a' = ";row_n(1)
print
mat print zzz
!****************************************************************
!combine = 'second column position' and 'thridth column position'
!****************************************************************
for i1=2 to col
call del_col_matrix(zzz(,),aaa(,),row_n(i1)-(i1-1),1)
mat print aaa
call resize_mat(zzz(,),aaa(,),nx,ny-(i1))
print " result :: delete column 'a' = ";row_n(i1)
print
mat print zzz
next i1
print "**************************************************************"
print " result = 'delete same : columns and rows ' from ' matrix a ' "
print "**************************************************************"
print
dim result(1,1)
mat redim result(size(zzz,2),size(zzz,1))
mat result = trn(zzz)
print "*****************************"
print " result ( again transpose) = "
print "*****************************"
mat print result
! data ( read - data statement)
data 1,2,3,4,5
data 6,7,8,9,10
data 11,12,13,14,15
data 16,17,18,19,20
data 21,22,23,24,25
end
sub resize_mat(zz(,),aa(,),nx,ny)
!*********************************************
! size reduction of the matrix (without zeros)
!*********************************************
mat redim zz(nx,ny)
mat zz = 0
for i=1 to nx
for j=1 to ny
zz(i,j) = aa(i,j)
next j
next i
end sub
sub del_col_matrix(a(,),aa(,),col_n,nc)
! fill (without size reduction)
for i=1 to size(a,1)
for j=1 to size(a,2)-1
if j < col_n then
aa(i,j) = a(i,j)
else
aa(i,j) = a(i,j+1)
end if
next j
next i
end sub
output view : deletecolumnsrowssplitscreen
given: three sides a,b,c ( triangle ?),type(obtuse,right,acute :triangle),option
---------------------------------------------------------------------------------------------------------------------
option nolet
!---------------------------------------
! find out what kind of triangle,option
! filename:TRIANGLETYPEoptionfinal.tru
!---------------------------------------
print "*******************************"
print " type : triangle(type),option "
print "*******************************"
print
print "**************************************"
print "[email protected],27/06/2016"
print "**************************************"
print
!*********************************
!** change by input statement ****
!*********************************
a=5 ! side 'a' of triangle
b=4 ! side 'b' of triangle
c=3 ! side 'c' of triangle
!*********************************
! tri = triangle
!*********************************
call print_info(a,b,c)
print
call check_tri_no_yes(a,b,c,flag)
print
print " flag = ";flag
print
call answer_no_yes(flag)
print
call type_tri_option(a,b,c)
end
sub print_info(a,b,c)
print " length of three sides (? triangle) "
print
print " length 'side a' ";a
print " length 'side b' ";b
print " length 'side c' ";c
end sub
sub check_tri_no_yes(a,b,c,flag)
if a<b+c and b<c+a and c<a+b then
flag =1
else
flag =0
end if
end sub
sub answer_no_yes(flag)
if flag = 0 then
print " No triangle forming witn sides a,b,c "
pause 5
stop
else
print " triangle forming with sides a,b,c "
end if
end sub
sub type_tri_option(a,b,c)
declare function perimeter
if a^2 > b^2+c^2 or b^2 > a^2+c^2 or c^2 > a^2+b^2 then
print " triangle :(type) is obtuse "
call option_type(a,b,c)
print " perimeter = ";perimeter(a,b,c)
call area_triangle(a,b,c,area)
print " area triangle = ";area
exit if
elseif a^2=b^2+c^2 or b^2=a^2+c^2 or c^2=a^2+b^2 then
! 2e) type right angle triangle
print " triangle:(type) right angle triangle "
call option_type(a,b,c)
print " perimeter = ";perimeter(a,b,c)
call area_triangle(a,b,c,area)
print " area triangle = ";area
exit if
else
! 3e) type acute angled triangle
print " triangle:(type) acute angle triangle "
call option_type(a,b,c)
print " perimeter = ";perimeter(a,b,c)
call area_triangle(a,b,c,area)
print " area triangle = ";area
exit if
end if
end sub
function perimeter(a,b,c)
perimeter = a+b+c
end function
sub area_triangle(a,b,c,area)
! perimeter = a+b+c
s=(a+b+c)/2
area = sqr(s*(s-a)*(s-b)*(s-c))
end sub
sub option_type(a,b,c)
if a=b and b=c then
print " option : triangle is equilateral"
elseif a=b or b=c or a=c then
print " option : triangle is isosceles "
else
print " option : triangle is scalene "
end if
end sub
output view : triangleyesnotypeoption
draw a trianglemesh (version 1.0)
-----------------------------------------------------
option nolet
!draw triangle mesh on the screen
!filename:DRAWTRIANGLEMESHFROMDATAfinalversion1.tru
dim xy(1,2),tri(1,3)
mat redim xy(5,2),tri(3,3)
! open #n:screen xl,xh,yl,yh
open #1:screen 0,0.5,0,1
open #2:screen .5,1,.5,1
window #1
print
print "***************************************"
print "* draw triangle mesh (version 1.0) *"
print "* [email protected],30/06/16*"
print "***************************************"
print
print "----------------------"
print " nodal values (x,y) : "
print "----------------------"
print
mat read xy
print " x"," y"
print
mat print xy
print "-------------------"
print " nodes triangles : "
print "-------------------"
print " local numbers "
print " 1"," 2"," 3"
mat read tri
print
mat print tri
window #2
ask pixels px,py
! aspect_ratio (px > py) => (px/py > 1)
asp_r = px/py
!****************
! *** x-range ***
x_low =-.5
x_high =7.2
!****************
! *** y-range ***
y_low =-.5
y_high =7.2
!****************
set window x_low*asp_r,x_high*asp_r,y_low,y_high
! draw axes
plot lines: x_low*asp_r,0;x_high*asp_r,0
plot lines: 0,y_low;0,y_high
plot text,at 0,0:"("& str$(0)&"," & str$(0) & ")"
! draw triangle :
plot lines: xy(1,1),xy(1,2);xy(2,1),xy(2,2)
plot lines: xy(2,1),xy(2,2);xy(3,1),xy(3,2)
plot lines: xy(1,1),xy(1,2);xy(3,1),xy(3,2)
!*****************************
! draw three 'triangles = tri'
!*****************************
for i=1 to size(tri,2)
call draw_tri(i,xy(,),tri(,))
next i
! draw node number
set color red
for i=1 to size(xy,1)
call node_nu(i,xy(,))
next i
! draw triangle number
dim xcyc(1,2)
nc=size(tri,2)
mat redim xcyc(nc,2)
set color blue
for i1=1 to nc
call tri_nu(i1,xy(,),tri(,),xcyc(,))
next i1
!************************
!* data - read statement*
!************************
!
!************************
! read xy 'nodal values'
!***********************
data 1,1,3,.5,6,1.5,3,2,7,2
!******************************
! read tri 'nodes:triangle',ccw
! ccw = counter clock wise
!******************************
data 1,2,3
data 1,3,4
data 3,5,4
end
sub tri_nu(i1,xy(,),tri(,),xcyc(,))
xcyc(i1,1)=(xy(tri(i1,1),1)+xy(tri(i1,2),1)+xy(tri(i1,3),1))/3
xcyc(i1,2)=(xy(tri(i1,1),2)+xy(tri(i1,2),2)+xy(tri(i1,3),2))/3
set color blue
plot text, at xcyc(i1,1),xcyc(i1,2):"("& str$(i1) & ")"
end sub
sub node_nu(i1,xy(,))
plot text, at xy(i1,1)+0.01,xy(i1,2):str$(i1)
end sub
sub draw_tri(i1,xy(,),tri(,))
! draw three lines forming a triangle
plot lines: xy(tri(i1,1),1),xy(tri(i1,1),2);xy(tri(i1,2),1),xy(tri(i1,2),2)
plot lines: xy(tri(i1,2),1),xy(tri(i1,2),2);xy(tri(i1,3),1),xy(tri(i1,3),2)
plot lines: xy(tri(i1,1),1),xy(tri(i1,1),2);xy(tri(i1,3),1),xy(tri(i1,3),2)
end sub
output view : drawtrianglemeshversion1
===================================================
=====================================================
Thiele-Geddes sumproducts
download code: Thielegeddessumproducts
output view : view Thielegeddessumproducts
option nolet
print"**************************************************************"
print" example | ( n=3 ) input: A0=1;A1=2;A2=3;A3=4 "
print"--------- "
print" called : Thiele-Geddes sumproducts "
print"**************************************************************"
print" sum = A0*A1*A2*A3+A1*A2*A3+A2*A3+A3+1=65 "
print"**************************************************************"
dim a(0 to 1)
dim b(0 to 1)
n=3 ! to change n you can make any computation (with this values above)
mat redim a(0 to n)
mat redim b(0 to n)
!**************************
! input values A0,A1,A2,A3*
!**************************
for i=0 to n
a(i)=i+1
next i
!-----------------------
! products:A0*A1*....*Ai
!-----------------------
prod=1 !init
for i=0 to n
prod=1
for j=i to n
prod=prod*a(j)
next j
b(i)=prod
next i
print "******************"
print "matrix of products"
print "******************"
mat print b
print"!---------------------"
print"! (sum of products) +1"
print"!---------------------"
som=1 !init
for i=0 to n
som=som+b(i)
next i
print "sum = ";som
end
graphics : sub routines (below) , general use when draw graphics on the screen.
sub draw_asp(xwin,ywin,asp)
set window -xwin,xwin,-ywin,ywin
! control aspect ratio for circle not ellips
ask pixels px,py
if px>py then
asp=px/py
else
asp=py/px
end if
end sub
SUB draw_axes(xmin,xmax,ymin,ymax,nt)
LET ntick = nt ! number of tick marks
! distance between tick marks on x-axis
LET dx = (xmax - xmin)/ntick
! distance between tick marks on y-axis
LET dy = (ymax - ymin)/ntick
! include margin
SET WINDOW xmin - dx,xmax + dx,ymin - dy,ymax + dy
LET x0 = max(0,xmin)
LET y0 = max(0,ymin)
IF ymin*ymax < 0 then
LET y0 = 0
ELSE
LET y0 = ymin
END IF
PLOT LINES: xmin,y0;xmax,y0 ! horizontal axis
PLOT LINES: x0,ymin;x0,ymax ! vertical axis
LET tx = 0.1*dy ! size of tick mark
LET ty = 0.1*dx
FOR itick = 0 to ntick
LET x = xmin + itick*dx
PLOT LINES: x,y0 - tx; x,y0 + tx ! draw ticks on x axis
LET y = ymin + itick*dy
PLOT LINES: x0 - ty,y; x0 + ty,y ! draw ticks on y axis
NEXT itick
END SUB
sub smallcir_dr(s1,xf,yf,scale,asp,z3$)
set color z3$
for theta=0 to 2*3.14 step 0.01
plot points:s1*cos(theta)+xf,s1*sin(theta)*asp+yf
next theta
end sub
================================================================
The use of the 'mat' command.
1e) sub
sub somprodmat(i,j,m(,),n(,),c(,),bg,be)
som=0
for k=bg to be
som = som + m(i,k)*n(k,j)
next k
c(i,j) = som
end sub
sub somprodvec(i,m(,),z(),cc(),bg,be)
som=0
for k=bg to be
som = som + m(i,k)*z(k)
next k
cc(i) = som
end sub
2e) 'see output view : below ,sub below'
option nolet
print "*****************************"
print " info ,'mat' command "
print " [email protected] "
print "*****************************"
dim a(2,1)
a(1,1)= 7
a(2,1)=8
dim b(1,2)
b(1,1)=5
b(1,2)=6
dim c(2,2)
c(1,1)=1
c(1,2)=2
c(2,1)=3
c(2,2)=4
dim e(2,2)
e(1,1)=1
e(1,2)=2
e(2,1)=3
e(2,2)=4
dim d(1,2)
print " given = matrix b,c "
print " b= "
mat print b
print " c= "
mat print c
!*************************************
! multiplication two matrices
! b(1,2)*c(2,2) => d(1,2)
! (1,2)*(2,2) contract to (1,2)
!*************************************
! 1e) (1,2)*(2,2) = (1,2)
!*************************************
call prod_mat(1,2,2,b(,),c(,),d(,))
print " matrix = d =>multiplication matrix 'b*c' "
print " dimension => dim d(1,2) "
mat print d
dim cmat(1,1)
!*************************************
! multiplicatication two matrices
! d(1,2)*a(2,1) = > number
!*************************************
! 2e) (1,2)*(2,1) = (1,1) = number
!*************************************
print " given = matrix a "
print " a= "
mat print a
print "matrix = cmat =>multiplication matrix 'd*a'"
call prod_mat(1,2,1,d(,),a(,),cmat(,))
print "result => matrix multiplication = 'b*c*a' "
print " cmat ' number ' = "
mat print cmat
kk=1
ll=-4
mat redim cmat(2,2)
i=2
j=2
print " sum of two matrix kk= ";kk;" ll= ";ll
print " given = matrix c,e "
print "c= "
mat print c
print "e= "
mat print e
call som_mat(kk,ll,i,j,c(,),e(,),cmat(,))
print " result cmat = ";kk;"*c+ ";ll;"*e"
mat print cmat
end
sub prod_mat(i,k,j,amat(,),bmat(,),cmat(,))
mat redim amat(i,k),bmat(k,j),cmat(i,j)
mat cmat = amat*bmat
end sub
sub som_mat(kk,ll,i,j,amat(,),bmat(,),cmat(,))
!use :kk,ll float
! i,j integer
mat redim amat(i,j),bmat(i,j),cmat(i,j)
mat amat =kk*amat
mat bmat =ll*bmat
mat cmat = amat+bmat
end sub
==================================================================
download code: basicmatrixconstruct
output view :viewbasicmatrixconstruct
special functions ( example : series : Bessel function of order n (integer) = Jn(t))
n = integer include 0,1,2,3,....
t = ?
coef = t^n/(2^n*n!)
=> Jn(t) = coef*(1-t^2/(2*(2*n+2))+t^4/(2*4*(2*n+2)*(2*n+4))-.........)
option nolet
! special functions(Legrendre,Bessel,gamma......"theoretical chemistry=schrodinger eq")
print "******************************************************"
print " chemical engineering "
print " Bessel n th order (integer) version 1.0 "
print " [email protected] , 24-11-2015 "
print "******************************************************"
!-------------------------------------------------------------
! ref : schaum's outline series
! laplace transforms !(time => freq domain)
! !( inverse laplace transform ' freq => time domain ')
! Murray R.SPIEGEL
! McGraw-Hill BOOK COMPANY
! including 450 solved problems ( page 7,edition 1965)
! general : every edition can be used
!------------------------------------------------------------
!************************************************************
! 1e)zero order besselfunction. = b0 = f_Bessel0
! 2e) n (integer order) besselfunction = bn = b_n
!************************************************************
! example =>
! serie infinite number of terms (4 terms here => q=4,t= ?)
!************************************************************
t=1
t2=t*t
t4=t2*t2
t6=t4*t2
b0 = 1-t2/4+t4/(2*2*4*4)-t6/(2*2*4*4*6*6)
print " b0 (bessel zero order) = ";b0
!*****************************************************
q=3 ! first term ( som = 1 ) init => q=4
!*****************************************************
! programming series :bessel n(integer) order =>
!*****************************************************
call f_Bessel0(t,q,som)
print " f_Bessel0 (bessel zero order )= ";som
p=0 ! integer
! without coef
p2 = 2*(2*p+2)
p4 = p2*4*(2*p+4)
p6 = p4*6*(2*p+6)
bn = 1-t2/p2+t4/p4-t6/p6
print " bn (integer) = ";bn
call f_Besseln(p,t,q,som)
print " bn (integer) = ";som
pp = 0
call faculteit(pp,pr)
print " pr = ";pr
call coefbesseln(t,pp,pr,coef_besseln)
b_n = som*coef_besseln
print " with coef besseln (bessel n(integer) order = ";b_n
end
sub f_Bessel0(t,q,som)
product=1
som =1 ! first term is 1 otherwise use 0
for m=1 to q
product = 1
for n=1 to m
product = product*(2*n)*(2*n)
next n
som = som + ((-1)^m*t^(2*m))/product
next m
end sub
sub f_Besseln(p,t,q,som)
! without coef
product=1
som =1 ! first term is 1 otherwise use 0
for m=1 to q
product = 1
for n=1 to m
product = product*(2*n)*(2*p+(2*n))
next n
som = som + ((-1)^m*t^(2*m))/product
next m
end sub
sub faculteit(pp,pr)
if pp = 0 then
pr= 1
else
pr = 1
for i=1 to pp
pr=pr*i
next i
end if
end sub
sub coefbesseln(t,pp,pr,coef_besseln)
coef_besseln = t^pp/(2^pp*pr)
end sub
=======================================================
def :kgv=kleinste gemeen veelvoud(= dutch) = lcd( least common divider)
ggd=grootste gemene deler (= dutch) = gcd( greatest common divider)
lemma : ggd(a,b)*kgv(a,b)=a*b
option nolet
declare function ggd,kgv
a=3
b=26
c=a*b
ant1=ggd(a,b)
print "grootste gemene deler(";a;",";b;") = ";ant1
ant2=kgv(ant1,a,b)
print "kleinste gemene veelvoud(";a;",";b;") = ";ant2
! floor are the same int
f=int(2.501)
print f
end
def kgv(d,a,b)
c=a*b
kgv=c/d
! ggd(a,b)*kgv(a,b)=a*b
end function
def ggd(d,e)
do
tmp = mod(d,e)
d = e
e = tmp
loop while e <> 0
ggd = d
end function
! gcd simple (greatest common divisor)
! find the common factors of two numbers
! number m,number n
! print "list: use matrix command"
! find: gcd(m,n)
option nolet
m = 63
n = 14
print "***********************************"
print " peter.vlasschaert ,26/07/2016 "
print " gcd:(g)reatest (c)ommon (d)ivisor "
print "***********************************"
print " m = number ";m
print " n = number ";n
dim fm(1),fn(1)
!fm=factors of m: append to list fm
n1=0
for i=1 to m
if mod(m,i)=0 then
n1=n1+1
mat redim fm(n1)
fm(n1)=i
end if
next i
print " factors of fm = "
mat print fm
!***********************************
! result fm : 1,3,7,9,21,63
!***********************************
!fn=factors of n: append to list fn
m1=0
for i=1 to n
if mod(n,i)=0 then
m1=m1+1
mat redim fn(m1)
fn(m1)=i
end if
next i
print " factors of fn = "
mat print fn
!***********************************
! result fm : 1,2,7,14
!***********************************
dim cf(1)
! cf=common factors of " fm 'and' fn "
! at least one common divisor '1'
l1=0
n_l1=size(fn)
m_l1=size(fm)
zn = min(m_l1,n_l1)
for i=1 to zn
if fm(i)=fn(i) then
l1=l1+1
mat redim cf(l1)
! cf(l1) = fn(i),choose this one is also possible
cf(l1)=fm(i)
end if
next i
print "list 'cf' of common factors of 'fm and fn' "
mat print cf
!*****************
! result cf:1,7
!*****************
!
! find highest common factor (last in the list)
!
n_cf = size(cf)
gcd = cf(n_cf)
print " gcd(";m;",";n;")= ";gcd
!********************
! result gcd(63,14)=7
!********************
end
recursion (Wikipedia),factorial,sinus
---------------------------------------------------------
see above: other way to get factorial definition.
' sub faculteit(pp,pr) '
!*****************************************************
! recursion : factorial function
! fac(n) = 1*2*3...(n-1)*n
! recursion : sinus function ? How see below
! sin(x) = x - x^3/fac(3)+x^5/fac(5)-x^7/fac(7)+...
! recursion : cosinus function ? How see below
! cos(x) = 1 - x^2/fac(2)+x^4/fac(4)-x^6/fac(6)+
!*****************************************************
print "*********************************************"
print " recursion : sin function ,cos function 1.1 "
print " [email protected] , 17/05/2016 "
print " ********************************************"
option nolet
declare function fac,a,b
print "fac(";0;")=";fac(0)
for k = 0 to 3
ii = 2*k+1
print "fac(";ii;")=";fac(ii)
next k
print " sin(pi/4) = ";sin(pi/4)
print " sin(3.14159/4) = "; sin(3.14159/4)
x = 3.14159/4
print " x = ";x
nu = 4 ! first term = x;second term = - x^3/fac(3),.....
print " number of terms for 'recursion' = ";nu
print "*************************************************"
print " sin(x) = x - x^3/fac(3)+x^5/fac(5)-x^7/fac(7) *"
print "*************************************************"
print " sin(pi/4) = ";x-x^3/fac(3)+x^5/fac(5)-x^7/fac(7)
print "********************"
print " recursion : sin(x) "
print "********************"
print "*************************************************"
print " t(k) = (-1)^(k-1)*x^(2*k-1)/fac(2*k-1) "
print " next term : k -> k+1 "
print " t(k+1) = (-1)^k*x^(2*k+1)/fac(2*k+1) "
print "*************************************************"
print " t(k+1) = a(k,x) * t(k) "
print " ************************************************"
print " a(k,x) = t(k+1)/t(k) = - x^2/((2*k)*(2*k+1)) "
print "*************************************************"
sum_sin = x ! first term 'sinus serie '
term = x ! first term
for k = 1 to nu-1 ! nu-1=k => t(k+1)=t(nu) terms
! recursion ' inside '
term = a(x,k)*term
sum_sin = sum_sin + term
next k
print "recursion : sin(";x;") = ";sum_sin
print " sin(";x;") = sqr(2)/2 = ";sqr(2)/2
for k = 1 to 3
ii = 2*k
print "fac(";ii;")=";fac(ii)
next k
print "***********************************************"
print " cos(x) = 1 - x^2/fac(2)+x^4/fac(4)-x^6/fac(6) "
print "***********************************************"
print " (sum_cos)^2+(sum_sin)^2 = 1 'identity' "
print "***********************************************"
sum_cos = sqr(1-(sum_sin)^2)
print " cos(";x;") = ";sum_cos
print " serie cos(";x;")= ";1 - x^2/fac(2)+x^4/fac(4)-x^6/fac(6)
sum_cos = 0 ! classic ' sum = 0 '
term = 1 ! first term
for k = 1 to nu-1 ! nu-1=k => t(k+1)=t(nu) terms
! recursion ' inside '
sum_cos = sum_cos + term
term = b(x,k)*term
next k
print "********************"
print " recursion : cos(x) "
print "********************"
print " ************************************************"
print " b(k,x) = t(k+1)/t(k) = - x^2/((2*k)*(2*k-1)) "
print "*************************************************"
print " recursion : cos(";x;") = ";sum_cos
print " cos(";x;") = sqr(2)/2 = ";sqr(2)/2
end
def a(x,k)
a = - x^2/((2*k)*(2*k+1))
end def
def b(x,k)
b = - x^2/((2*k)*(2*k-1))
end def
def fac(i)
!***********************************************************
! recursion ' inside '
! https://en.wikipedia.org/wiki/Recursion_(computer_science)
!***********************************************************
if i = 0 then
fac = 1
else
fac = i*fac(i-1)
end if
end def
output view : view recursionsincosfunc1point1
TableForm : mat (usefull for finite elements )
-----------------------------------------------------------------------
! put smallmatrix in bigmatrix (for every element)*
! very handy tool for 2d matrix (table form) *
!**************************************************
! [email protected] , 30/05/2016 *
!**************************************************
print "**********************"
print "* tableform : matrix *"
print "**********************"
print " n = k1*(i-1) + j "
print "*************************"
print " example: k1 =5 "
print "*************************"
print " i,j =counter big matrix "
print "*************************"
!*******************************************************
! why : you put a other matrix in column 8 and then make
! the sum of every element from column 7 and column 8
! and put the sum in column 9 (you can do it )
!*******************************************************
option nolet
pos1 = 3
pos2 = 1
pos3 = 4
print " position => "
print "position(1)= "; pos1
print "position(2)= "; pos2
print "position(3)= "; pos3
print " matrix a => "
dim a(3,3)
!values matrix a ( 'read - data' statement)
for i=1 to 3
for j=1 to 3
read a(i,j)
next j
next i
mat print a
! dim ' big matrix '
k1 = 5
k2 = 5
k3 = k1*k2
dim k(1,1)
mat redim k(k1,k2)
mat k = zero
dim table(1,1)
mat redim table(k3,9)
! column 1
for i=1 to k3
table(i,1) = i
next i
! column 3,4
for i=0 to k1-1
for j=1 to k1
table(j+i*k1,3)=i+1
table(j+i*k1,4)=j
next j
next i
! algo
for i=1 to k3
if i=pos1 or i=pos2 or i=pos3 then !test: row ( position mat :a)
for j=1 to k3
if j=pos1 or j=pos2 or j=pos3 then !test: column ( position mat :a)
table(k1*(i-1)+j,5) = i ! column 5
table(k1*(i-1)+j,6) = j ! column 6
table(k1*(i-1)+j,2) = k1*(i-1)+j ! column 2
!*******************************************************************************
!use: example simplex 2d in a bigger matrix (finite element method:flow problem)
! : triangles ( simplex 2d),three points.
!*******************************************************************************
s1=a(1,1)
!*(pos1,pos1)*
k(pos1,pos1)=s1
table(k1*(pos1-1)+pos1,7) = s1 !column 7
s2=a(1,2)
!*(pos1,pos2)*
k(pos1,pos2)=s2
table(k1*(pos1-1)+pos2,7) = s2 !column 7
s3=a(1,3)
!*(pos1,pos3)*
k(pos1,pos3)=s3
table(k1*(pos1-1)+pos3,7) = s3 !column 7
s4=a(2,1)
!*(pos2,pos1)*
k(pos2,pos1)=s4
table(k1*(pos2-1)+pos1,7) = s4 !column 7
s5=a(2,2)
!*(pos2,pos2)*
k(pos2,pos2)=s5
table(k1*(pos2-1)+pos2,7) = s5 !column 7
s6=a(2,3)
!*(pos2,pos3)*
k(pos1,pos3)=s6
table(k1*(pos2-1)+pos3,7) = s6 !column 7
s7=a(3,1)
!*(pos3,pos1)*
k(pos3,pos1)=s7
table(k1*(pos3-1)+pos1,7) = s7 !column 7
s8=a(3,2)
!*(pos3,pos2)*
k(pos3,pos2)=s8
table(k1*(pos3-1)+pos2,7) = s8 !column 7
s9=a(3,3)
!*(pos3,pos3)*
k(pos3,pos3)=s9
table(k1*(pos3-1)+pos3,7) = s9 !column 7
end if
next j
end if
next i
print " table => "
print "number" ,"number"," i"," j","position i","position j","value:matrix a"
print 1,2,3,4,5,6,7,8,9,"column"
mat print table
data .3,.14,99
data .5,.8,.6
data .1,.11,.3
end
output view : tableformmat
TableForm : mat (relation : between 'index and number' )
-----------------------------------------------------------------------------------------
! find relation between nu=number and index i,j*
!***********************************************
! [email protected], 30/05/2016
!***********************************************
option nolet
print "******************************************************"
print " tableform : relation between ' index i,j and number' "
print "******************************************************"
maxj = 3 ! max length of the block
dim table(1,1)
k=maxj*maxj
mat redim table(k,4) ! make a table
mat table = zero ! make all entries of mat table zero
for i=1 to k
table(i,1) = i ! column '1' number 'nu'
next i
for i=0 to maxj-1
for j=1 to maxj
table(j+i*maxj,2)=i+1 ! column '2' index 'i'
table(j+i*maxj,3)=j ! column '3' index 'j'
next j
next i
print "display table : lengthblock = maxj = ";maxj
print " number = nu ","index 'i'","index 'j'"
mat print table
print "*************************************"
print "*************************************"
print "* nu = maxj*(i-1)+j Formula (1)"
print "*************************************"
print "*************************************"
ii = 3
print " given = i => ";ii
jj = 2
print " given = j => ";jj
nu1 = maxj*(ii-1)+jj
print " use Formula (1) 'see table to veryfy' "
print " find = nu => ";nu1
print "**************************************"
print "**************************************"
print " find inversion for 'Formula (1)' "
print "**************************************"
print " algo : 'Formula (2)' "
print " index 'j' => "
print " j = nu mod maxj "
print " if j = 0 then j = maxj "
print " if j <> 0 then j = nu mod maxj "
print " index 'i' => "
print " i = (nu-j)/maxj + 1 "
print "**************************************"
print "**************************************"
nu2 = 9
print " given = nu => ";nu2
!******************************
! implementation of Formula (2)
!******************************
j = mod(nu2,maxj)
print " result j 'check' j=0 or j<>0 => ";j
if j = 0 then
j = maxj
else
j=mod(nu2,maxj)
end if
print " check table to veryfy => "
print " find = j => ";j
i = (nu2-j)/maxj + 1
print " find = i => ";i
end
output view : tablerelationnumberindex
How to implement BC ( = Dirichlet conditions ) on a Stifness Matrix.
--------------------------------------------------------------------------------------------------------
! how to implement boundary conditions (stifness matrix = 'k')
! => constant values ' start square system '
! first we have system of equations ( Ax=b)
! => finite elements : use A = 'k' (stifness matrix )
! : use b = 'p' (load vector)
! : use x = 'phi' (physical quantity )
! stifness matrix ' read - data '
!************************************************************
option nolet
! filename = boundaryconditiondirichletversion1.tru
print "*********************************************"
print "* how to implement BC (Dirichlet conditions)*"
print "* on a system of equations version 1.0 *"
print "*********************************************"
print "* [email protected],06/06/2016 *"
print "*********************************************"
read nu ! number of equations
dim k(1,1)
mat redim k(nu,nu)
for i=1 to nu
for j=1 to nu
read k(i,j)
next j
next i
print " stifness matrix => read - data "
mat print k
dim p(1,1)
mat redim p(nu,1)
print " load vector => read - data "
for i=1 to nu
read p(i,1)
next i
mat print p
read nc ! number of constraint
dim phi(1)
mat redim phi(nu) ! phi solution vector ( with constraint)
dim tt(1),tt_val(1)
mat redim tt(nc),tt_val(nc)
for i = 1 to nc
read tt(i)
next i
for i = 1 to nc
read tt_val(i)
next i
!************************************
!*constraint can put solution vector*
!************************************
for i=1 to nc
phi(tt(i))=tt_val(i)
next i
print "-------------------------"
print " before solve the system "
print "-------------------------"
mat print phi
! reduce system : ne = nu - nc
ne = nu-nc ! constraint system
dim k_co(1,1)
mat redim k_co(ne,nu)
dim un_tt(1)
mat redim un_tt(nu)
for i=1 to nu
for ii=1 to ne
!**** tt = {1,4,5}********
if tt(ii) = i then
un_tt(i) = tt(ii)
else
exit if
end if
next ii
next i
!******************************
!******************************
! construct ' data-file '
!**** un_tt = {1,0,0,4,5,0}****
!******************************
!******************************
!*************************************
! constraint stifness matrix = k_co
!*************************************
mat print un_tt
ii=0
for i=1 to nu
if un_tt(i) = 0 then ! we need the zero entries
ii = ii+1 ! values must be put in different rows
for j=1 to nu ! ( all values for the columns)
! read dimension (6 by 6) to a (ne by 6)
k_co(ii,j)=k(i,j)
next j
end if
next i
print " constrained stifness matrix => "
mat print k_co
!*************************************
! constraint load vector = p_co
!*************************************
dim p_co(1,1)
mat redim p_co(ne,1)
ii=0
for i=1 to nu
if un_tt(i) = 0 then ! we need the zero entries
ii = ii+1 ! values must be put in different rows
for j=1 to 1 !( j=1 for load-vector)
! read dimension (6 by 6) to a (ne by 6)
p_co(ii,j)=p(i,j)
next j
end if
next i
print " constrained load vector => "
mat print p_co
!*************************************************
! system to solve = 'k_sys '
!*************************************************
dim k_sys(1,1)
mat redim k_sys(ne,nc)
s=0
for i=1 to nu
if un_tt(i) = 0 then
s=s+1
for ii=1 to ne
k_sys(ii,s)=k_co(ii,i)
! k_sys(ii,2)=k_co(ii,3)
! k_sys(ii,3)=k_co(ii,6)
next ii
else
exit if
end if
next i
print " only system to solve : after put in the boundary conditions "
print "k_sys => "
mat print k_sys
!*************************************************
! known values,put opposit (-) sign to load vector
!*************************************************
dim p_remove(1,1)
mat redim p_remove(ne,1)
for ii=1 to nc
for i=1 to ne
p_remove(i,1)=k_co(i,tt(ii))
p_remove(i,1)= - phi(tt(ii))*p_remove(i,1)
next i
! mat print p_remove
next ii
! make a zero vector ' p_zero ' contain = sum of vectors
dim p_zero(1,1)
mat redim p_zero(ne,1)
!mat p_zero = 0
!****************************
! we can use also direct p_co
!****************************
mat p_zero = p_co
for ii=1 to nc
for i=1 to ne
p_remove(i,1)=k_co(i,tt(ii))
p_remove(i,1)= - phi(tt(ii))*p_remove(i,1)
next i
mat p_zero = p_zero + p_remove
next ii
!**********************************
! result pp of system we need solve
!**********************************
mat p_co = p_zero
print " load vector for 'k_sys' "
mat print p_co
!solve ( example 'data - file ' phi(2),phi(3),phi(6)'
dim k_inv(1,1)
mat redim k_inv(ne,nu-nc)
mat k_inv = inv(k_sys)
! mat print k_inv
dim result(1,1)
mat redim result(ne,1)
mat result = k_inv*p_co
! example :given 'phi(1),phi(4),phi(5) '
print " solve k_sys => "
mat print result
print "**************************************"
print " solution vector ph(1) to ph(";nu;") "
print " 'k'=stifness matrix 'p'=load vector "
print "**************************************"
ii = 0
for i=1 to nu
if un_tt(i) = 0 then
ii = ii+1
phi(i)=result(ii,1)
else
exit if
! print ii,i,phi(i)
end if
next i
mat print phi
data 6 ! nu=number of equations
! data : stifness matrix ( example 6 by 6)
data 1000,10,20,30,40,50
data 10,2000,21,31,41,51
data 20,21,3000,32,42,52
data 30,31,32,4000,43,53
data 40,41,42,43,5000,54
data 50,51,52,53,54,6000
! load vector
data 100,110,120,130,140,150
! constraint
data 3 ! nc = number of constraint
data 1,4,5 ! matrix tt ( i from phi )
data 5,-7,0 ! matrix tt_val ( ph(i) = values )
end
output view : boundaryconditiondirichletversion1
Parse a string : a$ ( with values inside and some delimiter = '|,;')
---------------------------------------------------------------------------------------------------
'SUB version ' = parse string a$
download code: DELIMITERSTRINGFINALSUBFINAL
output view : parsestringsubfinal
option nolet
!***************************************
! [email protected],08/06/2016
!***************************************
! parse a string : a$ !
!***************************************
print "*************************"
print "*** parse a string a$ ***"
print "*************************"
print "****************************************"
print "*[email protected],08/06/2016*"
print "****************************************"
a$ = "1.54,1.2,3.56"
print " string a$ = "
print a$
dim tab$(3)
! selection
tab$(1) = "|"
tab$(2) = ","
tab$(3) = ";"
print " find what type delimiter used : ? ' |,; ' "
for j=1 to size(tab$)
p1=pos(a$,tab$(j))
if p1 <> 0 then
pa = j
end if
next j
print pa ,tab$(pa)
print " find position first delimiter in a$ "
for i=1 to len(a$)
if a$[i:i]=tab$(pa) then
pb = i
exit for ! only find one : the first
end if
next i
print pb,a$[pb:pb]
print " find all positions for delimiter in a$"
dim pc(1)
ii = 0
! len : used for the length of the string
for i=1 to len(a$)
if a$[i:i]=tab$(pa) then
ii = ii+1
mat redim pc(ii)
pc(ii) = i
! exit for
end if
next i
nu = size(pc,1)
print " number of delimiters = ";nu
print " positions of all delimiters = "
mat print pc
! search for numbers = nu+1,example three strings
print " strings = "
print "first string = "; a$[1:pc(1)-1]
print "second string = "; a$[pc(1)+1:pc(2)-1]
print "thirdth string = "; a$[pc(2)+1:len(a$)]
print " convert : strings to values = "
print "first number = "; val(a$[1:pc(1)-1])
print "second number = "; val( a$[pc(1)+1:pc(2)-1])
print "thirdth number = "; val( a$[pc(nu)+1:len(a$)])
print "***************************************************"
print "*make sum of ";nu+1;" numbers from the string a$ =*"
print "***************************************************"
! init sum
s=0
! first value
s=s+val(a$[1:pc(1)-1])
! between values
for i=1 to nu-1
s=s + val( a$[pc(i)+1:pc(i+1)-1])
next i
! last value
s=s+val( a$[pc(nu)+1:len(a$)])
print " the sum of three numbers from string = "
print " sum = ";s
end
output view : stringtonumbers
( read -write ) file: use principal 'string parse'
---------------------------------------------------------------------------------------------------
option nolet
!**************************
! data file : testexcel.dat
!**************************
! make this file :
!**************************
! 1,2,3,4
! 5,7,8,9
! 10,11,12,13
!**************************
print "*************************************************"
print " put in path : testexcel.dat "
print "example: C:\Prog\TrueBASIC GOLD v6\testexcel.dat "
print "*************************************************"
print " use : for finite elements (data-file) => "
filename$ = "C:\Prog\TrueBASIC GOLD v6\testexcel.dat"
! input prompt "name of file for data : ":filename$
dim a$(1),b$(1)
call read_file(filename$,a$())
call write_file(filename$,a$())
mat print a$
!************************************************
! use : above principal 'parse a string'
! example : three strings to parse (a$(i),i=1..3)
!************************************************
print "*********************"
print " content of the file "
print "*********************"
print " filename = ";filename$
for i=1 to size(a$)
print i;" line row = ";a$(i)
next i
!******************************
! output :
! 1 line row =1,2,3,4
! 2 line row =5,6,7,8
! 3 line row =9,10,11,12
!******************************
end
sub read_file(filename$,array$())
open #4 :name filename$,create newold,org text,access outin
i=0
do while more #4
i=i+1
mat redim array$(i)
line input #4:array$(i)
loop
close #4
end sub
sub write_file(filename$,array$())
open #4 :name filename$,create newold,org text,access outin
erase #4
set #4:margin maxnum
total = Ubound(array$)
for i=1 to total
print #4:array$(i)
next i
close #4
end sub
( write ) file: Use :Ax=b 'matrix and vector' to file
---------------------------------------------------------------------------
see :advanced programma's 67)
option nolet
dim a(1,1),b(1)
!*******************************************
! => Ax = b (system of equations)
!*******************************************
! filename :matrixvectorwritedatafinal.tru
print "************************************"
print " write data : vector = b,matrix = A "
print "************************************"
print " Use : export to other software : "
print " Tecplot ,Matlab,Scilab,Maxima ... "
print "************************************"
print " Use : Ax=b ,concept "
print "************************************"
print "*************************************"
print " [email protected],11/06/16"
print "*************************************"
mat redim a(5,5),b(5)
! read matrix A( read - data)
for i=1 to size(a,1)
for j=1 to size(a,2)
read a(i,j)
next j
next i
print " Show matrix A = "
mat print A
print "****************************************"
print " write to : "
print "****************************************"
print " ***************************************"
print "* C:\Prog\TrueBASIC GOLD v6\mat_a.dat *"
print "****************************************"
filename_mat$ = "C:\Prog\TrueBASIC GOLD v6\mat_a.dat"
call write_matrix(filename_mat$,a(,))
! read vector b (read - data)
print " Show vector b = "
for i=1 to Ubound(b)
read b(i)
next i
mat print b
print "****************************************"
print " write to : "
print "****************************************"
print " ***************************************"
print "* C:\Prog\TrueBASIC GOLD v6\vec_b.dat *"
print "****************************************"
filename_vec$ = "C:\Prog\TrueBASIC GOLD v6\vec_b.dat"
call write_vector(filename_vec$,b())
! data Matrix A (read-data)
data 5,2,-4,3,2
data 8,-6,5,1,6
data 4,3,-9,12.5,6
data 8,6,-5,2,8
data 2,3,4,-5,2
! data Vector b (read-data)
data -2,4,3,1,7
end
! lib write data to file(test windows 10)
Sub write_vector(filename$,data())
open #1: name filename$,create newold,org text,access outin
erase #1
set #1:Margin maxnum
mat print #1:data
close #1
end sub
Sub write_matrix(filename$,data(,))
! write matrix to file
open #1: name filename$,create newold,org text,access outin
erase #1
set #1:Margin maxnum
Mat data = TRN(data)
Mat print #1:data
close #1
end sub
( delete ) matrix: delete column and row from matrix a
--------------------------------------------------------------------------------------
option nolet
!****************************************
!****************************************
! delete columns and rows (version 1.0)
! [email protected] ,13/06/2016
!****************************************
!****************************************
! use : together normal and first alone
!****************************************
!****************************************
! select column => col = ? (first )
! select row => row = ? (second)
!****************************************
!****************************************
! deletecolumnandrowmatrixafinal.tru
dim a(1,1)
mat redim a(5,4)
mat read a
! dimension (number rows,number columns)
n = size(a,1)
m = size(a,2)
print "************"
print " matrix A = "
print "************"
mat print a
print "********************************"
print "* DELETE COLUMN AND ROW FROM A *"
print "********************************"
col = 2
del_c = col ! delete column 2
print "***********************************"
print "*operation:delete ";del_c;".column*"
print "***********************************"
dim v(1,1)
mat redim v(n,m-1)
print "**********************************"
print "* after :operation :col del *"
print "**********************************"
call delete_column(n,m,col,a(,),v(,))
mat print v
! v :operation delete row
! v -> w matrix
!**********************************
! ref : matrix a
!**********************************
! n = number rows matrix 'v'
! m-1 = number columns matrix 'v'
! result : delete row = row
! n-1 = number rows matrix 'w'
! m-1 = number columns matrix 'w'
!**********************************
row = 2
del_r =row ! delete row 2
print "***********************************"
print "* operation:delete";del_r;".row *"
print "***********************************"
dim w(1,1)
mat redim w(n-1,m-1)
print "**********************************"
print "* after :operation :row del *"
print "**********************************"
call delete_row(n,m,row,v(,),w(,))
mat print w
! data - mat read (statement)
data 1,2,3,4,5
data 6,7,8,9,10
data 11,12,13,14,15
data 16,17,18,19,20
end
!*******************************************
! sub :first operation then second operation
! a->v then v->w => result : a->w
!*******************************************
sub delete_row(n,m,row,v(,),w(,))
!**********************************
! second: operation 'v->w'
!**********************************
! n = number rows matrix 'v'
! m-1 = number columns matrix 'v'
! result : operation (delete row)
! n-1 = number rows matrix 'w'
! m-1 = number columns matrix 'w'
!**********************************
del_r = row
for j=1 to m-1
for i=1 to n-1
if i >= del_r then
w(i,j) = v(i+1,j)
else
w(i,j) = v(i,j)
end if
next i
next j
end sub
sub delete_column(n,m,col,a(,),v(,))
!**********************************
! first: operation ' a -> v '
!**********************************
! n= number rows matrix 'a'
! m= number columns matrix 'a'
! result : operation (delete column)
! n = number rows matrix 'v'
! m-1 = number columns matrix 'v'
!**********************************
del_c = col ! delete column = col
for i=1 to n
for j=1 to m-1
if j >= del_c then
v(i,j) = a(i,j+1)
else
v(i,j) = a(i,j)
end if
next j
next i
end sub
output view : matrixoperationdel col and row
( delete ) matrix: delete more columns and rows from matrix a
-------------------------------------------------------------------------------------------------
option nolet
!*********************************************
!*********************************************
! delete more columns and rows (version 1.0)
! [email protected] ,14/06/2016
!*********************************************
! moretimesDELETECOLUMNANDROWsimple45forfinal.tru
!*********************************************
!
!*********************************************
!*********************************************
dim a(1,1),col(1),row(1)
mat redim a(4,5),col(3),row(3)
!***************************************
! read - data statement
!***************************************
mat read col
mat read row
mat read a
!***************************************
!***************************************
! dimension (number rows,number columns)
n = size(a,1)
m = size(a,2)
print "************************"
print "* matrix A = *"
print "************************"
print "*****************************************"
print " operation mat : delete columns and rows "
print "*****************************************"
mat print a
dim v(1,1),w(1,1)
!*****************************
! first step : one column,row
! first run =>
!*****************************
print " first run => "
mat redim v(n,m-1)
call delete_column(n,m,col(1),a(,),v(,))
print " delete column = ";col(1)
mat print v
mat redim w(n-1,m-1)
call delete_row(n,m,row(1),v(,),w(,))
print " delete row = ";row(1)
mat print w
!***********************
! second thridth combine
!************************
for ii = 2 to size(col)
print ii;"(e the : run "
mat redim v(n-(ii-1),m-ii)
call delete_column(n-(ii-1),m-(ii-1),col(ii),w(,),v(,))
print " delete column = ";col(ii)
mat print v
mat redim w(n-ii,m-ii)
call delete_row(n-(ii-1),m-(ii-1),row(ii),v(,),w(,))
print " delete row = ";row(ii)
mat print w
next ii
!******************************
!****** read - data ***********
!******************************
! data - mat read (statement)*
!******************************
!******************************
!******************************
!column : col
data 2,3,2
! row : row
data 2,3,2
! mat a
data 1,2,3,4,5
data 6,7,8,9,10
data 11,12,13,14,15
data 16,17,18,19,20
end
!*******************************************
! sub :first operation then second operation
! a->v then v->w => result : a->w
!*******************************************
sub delete_row(n,m,row,v(,),w(,))
!**********************************
! second: operation 'v->w'
!**********************************
! n = number rows matrix 'v'
! m-1 = number columns matrix 'v'
! result : operation (delete row)
! n-1 = number rows matrix 'w'
! m-1 = number columns matrix 'w'
!**********************************
del_r = row
for j=1 to m-1
for i=1 to n-1
if i >= del_r then
w(i,j) = v(i+1,j)
else
w(i,j) = v(i,j)
end if
next i
next j
end sub
sub delete_column(n,m,col,a(,),v(,))
!**********************************
! first: operation ' a -> v '
!**********************************
! n= number rows matrix 'a'
! m= number columns matrix 'a'
! result : operation (delete column)
! n = number rows matrix 'v'
! m-1 = number columns matrix 'v'
!**********************************
del_c = col ! delete column = col
for i=1 to n
for j=1 to m-1
if j >= del_c then
v(i,j) = a(i,j+1)
else
v(i,j) = a(i,j)
end if
next j
next i
end sub
output view :moredeletecolrowsforversion1
( delete ) vector: delete elements from vector (sub version)
--------------------------------------------------------------------------------------------
option nolet
!***********************************************
!***********************************************
! delete : elements from vector (version 1.0)
! [email protected] , 15/06/2016
!***********************************************
!filename : MORETIMESDELETEvectorfinalforversionfinal
dim a(1),aa(1)
mat redim a(13),aa(13)
mat read a
print "********************************************"
print "* initial ' vector a '(operations,delete) = "
print "********************************************"
print " size( vector 'a') = ";size(a)
mat print a
dim el(1)
nc = 3 ! number constraint
mat redim el(nc)
mat read el
print "******************************"
print " constraint : delete elements "
print "******************************"
for kf = 1 to nc
call del_el_vec(kf,el(),a())
print " delete element from vector 'a' = ";el(kf)
mat print a
next kf
print "*****************************************"
print " result = delete elements from vector 'a'"
print "*****************************************"
print " size(vector 'a') = ";size(a)-nc
mat redim a(size(a)-nc) ! remove last three positions
print " new vector 'a' = "
mat print a
!******************************
!****** read - data ***********
!******************************
! data - mat read (statement)*
!******************************
!******************************
!******************************
!
!******************************
!* read vector *
!******************************
data 2,0,0,0,4,0,0,0,2,0,0,0,0
! read : delete el
data 4,8,13
end
sub del_el_vec(kf,el(),a())
for i=1 to el(kf)-2
a(i) = a(i)
next i
for i=el(kf)-1 to size(a)-1
a(i) = a(i+1)
next i
end sub
output view : deleteelementsfromvectorouput
( write ) vector,mat: with delimiter = "," ";" (sub version)
--------------------------------------------------------------------------------------------
option nolet
! write vector to file seperated by 'point comma' = ;
print "**************************************************"
print "**************************************************"
print "*write : matrix and vector to file with delimiter*"
print "**************************************************"
print "* [email protected] , 18/06/2016 *"
print "**************************************************"
print "**************************************************"
dim aa(1)
mat redim aa(4)
mat read aa
delimiter$ = ";"
a$ = "C:\Prog\TrueBASIC GOLD v6\testaa.dat"
print "*************************"
print "*write vector:( to file)*"
print "*************************"
print " write to file (string) = ";a$
print " delimiter (string) = ";delimiter$
print " vector write to file (screen) = "
mat print aa
call write_vec_delimiter(a$,delimiter$,aa())
print " ************************************"
print " testaa.dat : "
print " "
print " 1;2;3;4 "
print "*************************************"
!
!*****************
! read-data vector
!*****************
! write matrix to file seperated by 'comma' = ,
data 1,2,3,4
dim bb(1,1)
mat redim bb(3,4)
mat read bb
delimiter$ = ","
b$ = "C:\Prog\TrueBASIC GOLD v6\testbb.dat"
print "*************************"
print "*write matrix:( to file)*"
print "*************************"
print " write to file (string) = ";b$
print " delimiter (string) = ";delimiter$
print " matrix write to file (screen) = "
!**************
!very important
!**************
mat bb = TRN(bb) ! transpose
!**************
mat print bb
call write_mat_delimiter(b$,delimiter$,bb(,))
print " ************************************"
print " testbb.dat : "
print " "
print " 1,2,3,4 "
print " 5,6,7,8 "
print " 9,10,11,12 "
print "*************************************"
!******************
! read-data matrix
!******************
data 1,2,3,4
data 5,6,7,8
data 9,10,11,12
end
!******************************************
! lib write : vec and matrix with delimiter
!******************************************
! 1e) write : vector with delimiter
sub write_vec_delimiter(a$,delimiter$,a())
open #1: name a$,create newold,org text,access outin
erase #1
set #1:Margin maxnum
for i=1 to Ubound(a)-1
print #1:a(i);delimiter$;
next i
for i = Ubound(a) to Ubound(a)
print #1:a(i);
next i
close #1
end sub
! 2e) write : matrix with delimiter
sub write_mat_delimiter(b$,delimiter$,bb(,))
open #2: name b$,create newold,org text,access outin
erase #2
set #2:Margin maxnum
n1 = size(bb,1)
n2 = size(bb,2)
!***********************************
! reverse direction of reading
! => We use transpose bb
! normal : first = i other for j
!***********************************
! last element : i control by i loop
!***********************************
for j=1 to n2
for i=1 to n1
if i<>n1 then
print #2:bb(i,j);delimiter$;
else
print #2:bb(i,j);
end if
next i
!*************************************************
! statement :
! very important otherwise :every element one line
!*************************************************
print #2:
next j
close #2
end sub
datafile vector :TESTAA.dat
datafile matrix :TESTBB.dat
output view : writevecandmatwithdelimiter
( delete ) matrix: delete same columns and rows from matrix a
=> split screen.
-------------------------------------------------------------------------------------------------
option nolet
! delete columns from matrix
dim a(1,1),aa(1,1)
mat redim a(5,5),aa(5,5)
!*****************************
! split screen two equal parts
!*****************************
! filename :
!DELETECOLUMNSMATMATRIX55finalfinalrowfinalsplitscreen.tru
open #1: screen 0,0.5,0,1
open #2: screen .5,1,0,1
!*****************************
! open left part of the screen
!*****************************
window #1
print "****************************************"
print " delete :columns and rows => matrix 'a' "
print " [email protected],25/06/2016 "
print "****************************************"
set zonewidth 20
!********************
! read data
!********************
mat read a
!*******************
! initial matrix a
!*******************
print " matrix a = "
mat print a
!************************
!************************
! column deleting matrix
!************************
dim col_n(1)
mat redim col_n(3)
col_n(1) = 2
col_n(2) = 4
col_n(3) = 5
col = size(col_n)
nx = size(a,1)
ny = size(a,2)
print "nx = ";nx
print "ny = ";ny
print " deleting positions of columns "
mat print col_n
!***********************
!***********************
!***********************
mat aa = zero
mat print aa
dim zz(1,1)
!****************************************************************
! 'first column position'
! delete column = 2 => col_n(1)-0
!****************************************************************
call del_col_matrix(a(,),aa(,),col_n(1)+0,1)
!mat print aa
call resize_mat(zz(,),aa(,),nx,ny)
print " result :: delete column 'a' = ";col_n(1)
mat print zz
!****************************************************************
!combine = 'second column position' and 'thridth column position'
!****************************************************************
for i1=2 to col
call del_col_matrix(zz(,),aa(,),col_n(i1)-(i1-1),1)
mat print aa
call resize_mat(zz(,),aa(,),nx,ny-i1)
print " result :: delete column 'a' = ";col_n(i1)
mat print zz
next i1
!mat print zz
dim row(1,1)
mat redim row(size(zz,2),size(zz,1))
! reason :same operations on the rows
mat row = trn(zz)
!******************************
! open right part of the screen
!******************************
window #2
set zonewidth 20
! mat 'a' equivalence with mat 'row'
print " matrix 'a' for rows (columm definition) "
print " mat row = transpose(zz) "
mat print row
dim aaa(1,1),zzz(1,1)
mat redim aaa(size(row,1),size(row,2))
nx=size(row,1)
ny=size(row,2)
print "nx = ";nx
print "ny = ";ny
print " deleting :same row positions "
dim row_n(1)
mat redim row_n(col)
mat row_n = col_n
mat print row_n
! a:part one
call del_col_matrix(row(,),aaa(,),row_n(1),1)
mat print aaa
! b:part one
call resize_mat(zzz(,),aaa(,),nx,ny)
print " result :: delete column 'a' = ";row_n(1)
mat print zzz
!****************************************************************
!combine = 'second column position' and 'thridth column position'
!****************************************************************
for i1=2 to col
call del_col_matrix(zzz(,),aaa(,),row_n(i1)-(i1-1),1)
mat print aaa
call resize_mat(zzz(,),aaa(,),nx,ny-(i1))
print " result :: delete column 'a' = ";row_n(i1)
mat print zzz
next i1
print "**************************************************************"
print " result = 'delete same : columns and rows ' from ' matrix a ' "
print "**************************************************************"
dim result(1,1)
mat redim result(size(zzz,2),size(zzz,1))
mat result = trn(zzz)
print "*****************************"
print " result ( again transpose) = "
print "*****************************"
mat print result
! data ( read - data statement)
data 1,2,3,4,5
data 6,7,8,9,10
data 11,12,13,14,15
data 16,17,18,19,20
data 21,22,23,24,25
end
sub resize_mat(zz(,),aa(,),nx,ny)
!*********************************************
! size reduction of the matrix (without zeros)
!*********************************************
mat redim zz(nx,ny)
mat zz = 0
for i=1 to nx
for j=1 to ny
zz(i,j) = aa(i,j)
next j
next i
end sub
sub del_col_matrix(a(,),aa(,),col_n,nc)
! fill (without size reduction)
for i=1 to size(a,1)
for j=1 to size(a,2)-1
if j < col_n then
aa(i,j) = a(i,j)
else
aa(i,j) = a(i,j+1)
end if
next j
next i
end sub
output view : deletecolumnsrowssplitscreen
given: three sides a,b,c ( triangle ?),type(obtuse,right,acute :triangle),option
---------------------------------------------------------------------------------------------------------------------
option nolet
!---------------------------------------
! find out what kind of triangle,option
! filename:TRIANGLETYPEoptionfinal.tru
!---------------------------------------
print "*******************************"
print " type : triangle(type),option "
print "*******************************"
print "**************************************"
print "[email protected],27/06/2016"
print "**************************************"
!*********************************
!** change by input statement ****
!*********************************
a=5 ! side 'a' of triangle
b=4 ! side 'b' of triangle
c=3 ! side 'c' of triangle
!*********************************
! tri = triangle
!*********************************
call print_info(a,b,c)
call check_tri_no_yes(a,b,c,flag)
print " flag = ";flag
call answer_no_yes(flag)
call type_tri_option(a,b,c)
end
sub print_info(a,b,c)
print " length of three sides (? triangle) "
print " length 'side a' ";a
print " length 'side b' ";b
print " length 'side c' ";c
end sub
sub check_tri_no_yes(a,b,c,flag)
if a<b+c and b<c+a and c<a+b then
flag =1
else
flag =0
end if
end sub
sub answer_no_yes(flag)
if flag = 0 then
print " No triangle forming witn sides a,b,c "
pause 5
stop
else
print " triangle forming with sides a,b,c "
end if
end sub
sub type_tri_option(a,b,c)
declare function perimeter
if a^2 > b^2+c^2 or b^2 > a^2+c^2 or c^2 > a^2+b^2 then
print " triangle :(type) is obtuse "
call option_type(a,b,c)
print " perimeter = ";perimeter(a,b,c)
call area_triangle(a,b,c,area)
print " area triangle = ";area
exit if
elseif a^2=b^2+c^2 or b^2=a^2+c^2 or c^2=a^2+b^2 then
! 2e) type right angle triangle
print " triangle:(type) right angle triangle "
call option_type(a,b,c)
print " perimeter = ";perimeter(a,b,c)
call area_triangle(a,b,c,area)
print " area triangle = ";area
exit if
else
! 3e) type acute angled triangle
print " triangle:(type) acute angle triangle "
call option_type(a,b,c)
print " perimeter = ";perimeter(a,b,c)
call area_triangle(a,b,c,area)
print " area triangle = ";area
exit if
end if
end sub
function perimeter(a,b,c)
perimeter = a+b+c
end function
sub area_triangle(a,b,c,area)
! perimeter = a+b+c
s=(a+b+c)/2
area = sqr(s*(s-a)*(s-b)*(s-c))
end sub
sub option_type(a,b,c)
if a=b and b=c then
print " option : triangle is equilateral"
elseif a=b or b=c or a=c then
print " option : triangle is isosceles "
else
print " option : triangle is scalene "
end if
end sub
output view : triangleyesnotypeoption
draw a trianglemesh (version 1.0)
-----------------------------------------------------
option nolet
!draw triangle mesh on the screen
!filename:DRAWTRIANGLEMESHFROMDATAfinalversion1.tru
dim xy(1,2),tri(1,3)
mat redim xy(5,2),tri(3,3)
! open #n:screen xl,xh,yl,yh
open #1:screen 0,0.5,0,1
open #2:screen .5,1,.5,1
window #1
print "***************************************"
print "* draw triangle mesh (version 1.0) *"
print "* [email protected],30/06/16*"
print "***************************************"
print "----------------------"
print " nodal values (x,y) : "
print "----------------------"
mat read xy
print " x"," y"
mat print xy
print "-------------------"
print " nodes triangles : "
print "-------------------"
print " local numbers "
print " 1"," 2"," 3"
mat read tri
mat print tri
window #2
ask pixels px,py
! aspect_ratio (px > py) => (px/py > 1)
asp_r = px/py
!****************
! *** x-range ***
x_low =-.5
x_high =7.2
!****************
! *** y-range ***
y_low =-.5
y_high =7.2
!****************
set window x_low*asp_r,x_high*asp_r,y_low,y_high
! draw axes
plot lines: x_low*asp_r,0;x_high*asp_r,0
plot lines: 0,y_low;0,y_high
plot text,at 0,0:"("& str$(0)&"," & str$(0) & ")"
! draw triangle :
plot lines: xy(1,1),xy(1,2);xy(2,1),xy(2,2)
plot lines: xy(2,1),xy(2,2);xy(3,1),xy(3,2)
plot lines: xy(1,1),xy(1,2);xy(3,1),xy(3,2)
!*****************************
! draw three 'triangles = tri'
!*****************************
for i=1 to size(tri,2)
call draw_tri(i,xy(,),tri(,))
next i
! draw node number
set color red
for i=1 to size(xy,1)
call node_nu(i,xy(,))
next i
! draw triangle number
dim xcyc(1,2)
nc=size(tri,2)
mat redim xcyc(nc,2)
set color blue
for i1=1 to nc
call tri_nu(i1,xy(,),tri(,),xcyc(,))
next i1
!************************
!* data - read statement*
!************************
!
!************************
! read xy 'nodal values'
!***********************
data 1,1,3,.5,6,1.5,3,2,7,2
!******************************
! read tri 'nodes:triangle',ccw
! ccw = counter clock wise
!******************************
data 1,2,3
data 1,3,4
data 3,5,4
end
sub tri_nu(i1,xy(,),tri(,),xcyc(,))
xcyc(i1,1)=(xy(tri(i1,1),1)+xy(tri(i1,2),1)+xy(tri(i1,3),1))/3
xcyc(i1,2)=(xy(tri(i1,1),2)+xy(tri(i1,2),2)+xy(tri(i1,3),2))/3
set color blue
plot text, at xcyc(i1,1),xcyc(i1,2):"("& str$(i1) & ")"
end sub
sub node_nu(i1,xy(,))
plot text, at xy(i1,1)+0.01,xy(i1,2):str$(i1)
end sub
sub draw_tri(i1,xy(,),tri(,))
! draw three lines forming a triangle
plot lines: xy(tri(i1,1),1),xy(tri(i1,1),2);xy(tri(i1,2),1),xy(tri(i1,2),2)
plot lines: xy(tri(i1,2),1),xy(tri(i1,2),2);xy(tri(i1,3),1),xy(tri(i1,3),2)
plot lines: xy(tri(i1,1),1),xy(tri(i1,1),2);xy(tri(i1,3),1),xy(tri(i1,3),2)
end sub
output view : drawtrianglemeshversion1
===================================================
=====================================================