Posts - MaplePrimes

Stiffness matrix implementation of a finite element...

by: Adjal yassine 25

November 07 2019

0 5

>

MD11 := subs(a[0] = 1, a[1] = 0, a[2] = 0, a[3] = 0, v[i]):

>

MS11 := subs(a[0] = 0, a[1] = 0, a[2] = 0, a[3] = 0, theta[i]):

>

MT11 := subs(a[0] = 1, a[1] = 0, a[2] = 0, a[3] = 0, v[j]):

>

MR11 := subs(a[0] = 1, a[1] = 0, a[2] = 0, a[3] = 0, theta[j]):

>

$KG := Matrix(4, 4, {(1, 1) = 12*EI/L^3, (1, 2) = 6*EI/L^2, (1, 3) = -12*EI/L^3, (1, 4) = 6*EI/L^2, (2, 1) = 6*EI/L^2, (2, 2) = 4*EI/L, (2, 3) = -6*EI/L^2, (2, 4) = 2*EI/L, (3, 1) = -12*EI/L^3, (3, 2) = -6*EI/L^2, (3, 3) = 12*EI/L^3, (3, 4) = -6*EI/L^2, (4, 1) = 6*EI/L^2, (4, 2) = 2*EI/L, (4, 3) = -6*EI/L^2, (4, 4) = 4*EI/L})$

(1)

>

Download stiffnesmatrice.mw

Stiffness matrix implementation of a finite element beam based on Euler-Bernoulli theory stiffnesmatrice.mw

Free vibration of a cracked composite beam

by: Adjal yassine 25 Product: Maple

November 06 2019

5 0

>	restart: with(LinearAlgebra):

>	# Motion equation ( Vibration of a cracked composite beam using general solution) Edited by Adjal Yassine #

>	####################################################################

Motion equation of flexural vibration in normalized form

>	EIW^(iv)-momega^2*W=0;

(1)

The general solution form of bending vibration equation

>	W1:=A[1]cosh(mux)+A[2]sinh(mux)+A[3]cos(mux)+A[4]sin(mux);

(2)

where

>	E:=2682e6;L:=0.18;h:=0.004;b:=0.02;rho:=2600;area=bh;m:=rhohb;II:=(hb^3)/12:

(3)

>	mu:=((momega^2L^4/EI)^(1/4)):

Expression of cross-sectional rotation , the bending moment shear force and torsional moment are given as follows respectively

>	theta1 := (1/L)(A[1]musinh(mux)+A[2]mucosh(mux)-A[3]musin(mux)+A[4]mucos(mu*x));

(A[1]*mu*sinh(mu*x)+A[2]*mu*cosh(mu*x)-A[3]*mu*sin(mu*x)+A[4]*mu*cos(mu*x))/L

(4)

>	M1:= (EI/L^2)(A[1]mu^2cosh(mux)+A[2]mu^2sinh(mux)-A[3]mu^2cos(mux)-A[4]mu^2sin(mu*x));

EI*(A[1]*mu^2*cosh(mu*x)+A[2]*mu^2*sinh(mu*x)-A[3]*mu^2*cos(mu*x)-A[4]*mu^2*sin(mu*x))/L^2

(5)

>	S1:= (-EI/L^3)(A[1]mu^3sinh(mux)+A[2]mu^3cosh(mux)+A[3]mu^3sin(mux)-A[4]mu^3cos(mu*x));

-EI*(A[1]*mu^3*sinh(mu*x)+A[2]*mu^3*cosh(mu*x)+A[3]*mu^3*sin(mu*x)-A[4]*mu^3*cos(mu*x))/L^3

(6)

>

>	W2:=A[5]cosh(mux)+A[6]sinh(mux)+A[7]cos(mux)+A[8]sin(mux);

(7)

>

>	theta2:= (1/L)(A[5]musinh(mux)+A[6]mucosh(mux)-A[7]musin(mux)+A[8]mucos(mu*x));

(A[5]*mu*sinh(mu*x)+A[6]*mu*cosh(mu*x)-A[7]*mu*sin(mu*x)+A[8]*mu*cos(mu*x))/L

(8)

>	M2:= (EI/L^2)(A[5]mu^2cosh(mux)+A[6]mu^2sinh(mux)-A[7]mu^2cos(mux)-A[8]mu^2sin(mu*x));

EI*(A[5]*mu^2*cosh(mu*x)+A[6]*mu^2*sinh(mu*x)-A[7]*mu^2*cos(mu*x)-A[8]*mu^2*sin(mu*x))/L^2

(9)

>	S2:= -(EI/L^3)(A[5]mu^3sinh(mux)+A[6]mu^3cosh(mux)+A[7]mu^3sin(mux)-A[8]mu^3cos(mu*x));

-EI*(A[5]*mu^3*sinh(mu*x)+A[6]*mu^3*cosh(mu*x)+A[7]*mu^3*sin(mu*x)-A[8]*mu^3*cos(mu*x))/L^3

(10)

>

The boundary conditions at fixed end W1(0)=Theta(0)=0

>	X1:=eval(subs(x=0,W1));

(11)

>	X2:=eval(subs(x=0,theta1));

(mu*A[2]+mu*A[4])/L

(12)

>	X2:=collect(X2,mu)*(L/mu);

(13)

>

The boundary condtions at free end M2(1)=S2(1)=0

>	X3:=eval(subs(x=1,M2));

EI*(A[5]*mu^2*cosh(mu)+A[6]*mu^2*sinh(mu)-A[7]*mu^2*cos(mu)-A[8]*mu^2*sin(mu))/L^2

(14)

>	X3:=collect(X3,mu)*(L^2/mu^2/EI);

(15)

>	X4:=eval(subs(x=1,S2));

-EI*(A[5]*mu^3*sinh(mu)+A[6]*mu^3*cosh(mu)+A[7]*mu^3*sin(mu)-A[8]*mu^3*cos(mu))/L^3

(16)

>	X4:=collect(X4,mu);

-EI*(cosh(mu)*A[6]+sinh(mu)*A[5]-cos(mu)*A[8]+sin(mu)*A[7])*mu^3/L^3

(17)

>	X4:=collect(X4,EI)*(L^3/mu^3/EI);

(18)

>

The additional boundary conditions at crack location

>	X5:=combine(M1-M2);

(EI*cosh(mu*x)*mu^2*A[1]-EI*cosh(mu*x)*mu^2*A[5]+EI*sinh(mu*x)*mu^2*A[2]-EI*sinh(mu*x)*mu^2*A[6]-EI*cos(mu*x)*mu^2*A[3]+EI*cos(mu*x)*mu^2*A[7]-EI*sin(mu*x)*mu^2*A[4]+EI*sin(mu*x)*mu^2*A[8])/L^2

(19)

>	X5:=collect(X5,mu);

(EI*cosh(mu*x)*A[1]-EI*cosh(mu*x)*A[5]+EI*sinh(mu*x)*A[2]-EI*sinh(mu*x)*A[6]-cos(mu*x)*EI*A[3]+A[7]*cos(mu*x)*EI-A[4]*sin(mu*x)*EI+A[8]*sin(mu*x)*EI)*mu^2/L^2

(20)

>	X5:=collect(X5,EI)*(L^2/mu^2/EI);

A[1]*cosh(mu*x)-A[5]*cosh(mu*x)+A[2]*sinh(mu*x)-A[6]*sinh(mu*x)-A[3]*cos(mu*x)+A[7]*cos(mu*x)-A[4]*sin(mu*x)+A[8]*sin(mu*x)

(21)

>	X6:=combine(S1-S2);

(-EI*cosh(mu*x)*mu^3*A[2]+EI*cosh(mu*x)*mu^3*A[6]-EI*sinh(mu*x)*mu^3*A[1]+EI*sinh(mu*x)*mu^3*A[5]+EI*cos(mu*x)*mu^3*A[4]-EI*cos(mu*x)*mu^3*A[8]-EI*sin(mu*x)*mu^3*A[3]+EI*sin(mu*x)*mu^3*A[7])/L^3

(22)

>	X6:=collect(X6,mu);

(-EI*cosh(mu*x)*A[2]+EI*cosh(mu*x)*A[6]-EI*sinh(mu*x)*A[1]+EI*A[5]*sinh(mu*x)+cos(mu*x)*A[4]*EI-cos(mu*x)*A[8]*EI-sin(mu*x)*EI*A[3]+sin(mu*x)*A[7]*EI)*mu^3/L^3

(23)

>	X6:=collect(X6,EI)*(L^3/mu^3/EI);

-cosh(mu*x)*A[2]+cosh(mu*x)*A[6]-sinh(mu*x)*A[1]+sinh(mu*x)*A[5]+cos(mu*x)*A[4]-cos(mu*x)*A[8]-sin(mu*x)*A[3]+sin(mu*x)*A[7]

(24)

>

>	X7:=combine(W2-(W1+c8*S1));

(EI*cosh(mu*x)*c8*mu^3*A[2]+EI*sinh(mu*x)*c8*mu^3*A[1]-EI*cos(mu*x)*c8*mu^3*A[4]+EI*sin(mu*x)*c8*mu^3*A[3]-cosh(mu*x)*A[1]*L^3+cosh(mu*x)*A[5]*L^3-sinh(mu*x)*A[2]*L^3+sinh(mu*x)*A[6]*L^3-cos(mu*x)*A[3]*L^3+cos(mu*x)*A[7]*L^3-sin(mu*x)*A[4]*L^3+sin(mu*x)*A[8]*L^3)/L^3

(25)

>	X8:=combine (theta2-(theta1+c44*M1));

(-EI*cosh(mu*x)*c44*mu^2*A[1]-EI*sinh(mu*x)*c44*mu^2*A[2]+EI*cos(mu*x)*c44*mu^2*A[3]+EI*sin(mu*x)*c44*mu^2*A[4]-L*cosh(mu*x)*mu*A[2]+L*cosh(mu*x)*mu*A[6]-L*sinh(mu*x)*mu*A[1]+L*sinh(mu*x)*mu*A[5]-L*cos(mu*x)*mu*A[4]+L*cos(mu*x)*mu*A[8]+L*sin(mu*x)*mu*A[3]-L*sin(mu*x)*mu*A[7])/L^2

(26)

>

The characteristic matrix function of frequency

>	FD8:=subs(A[1]=1,A[3]=0,X1);FD12:=subs(A[1]=0,A[3]=0,X1);FD13:=subs(A[1]=0,A[3]=1,X1);FD14:=subs(A[1]=0,A[3]=0,X1);FD15:=subs(A[1]=0,A[3]=0,X1);FD16:=subs(A[1]=0,A[3]=0,X1);FD17:=subs(A[1]=0,A[3]=0,X1);FD18:=subs(A[1]=0,A[3]=0,X1);

(27)

>	FD21:=subs(A[2]=0,A[4]=0,X2);FD22:=subs(A[2]=1,A[4]=0,X2);FD23:=subs(A[2]=0,A[4]=0,X2);FD24:=subs(A[2]=0,A[4]=1,X2);FD25:=subs(A[2]=0,A[4]=0,X2);FD26:=subs(A[2]=0,A[4]=0,X2);FD27:=subs(A[2]=0,A[4]=0,X2);FD28:=subs(A[2]=0,A[4]=0,X2);

(28)

>

FD31:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD32:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD33:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD34:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD35:=subs(A[5]=1,A[6]=0,A[7]=0,A[8]=0,X3);;FD36:=subs(A[5]=0,A[6]=1,A[7]=0,A[8]=0,X3);FD37:=subs(A[5]=0,A[6]=0,A[7]=1,A[8]=0,X3);FD38:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=1,X3);

(29)

>

FD41:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD42:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD43:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD44:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD45:=subs(A[5]=1,A[6]=0,A[7]=0,A[8]=0,X4);FD46:=subs(A[5]=0,A[6]=1,A[7]=0,A[8]=0,X4);FD47:=subs(A[5]=0,A[6]=0,A[7]=1,A[8]=0,X4);FD48:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=1,X4);

(30)

>

FD51:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD52:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD53:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD54:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD55:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X5);FD56:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X5);FD57:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X5);FD58:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X5);

(31)

>

FD61:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD62:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD63:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD64:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD65:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X6);FD66:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X6);FD67:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X6);FD68:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X6);

(32)

>

FD71:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD72:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD73:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD74:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD75:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X7);FD76:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X7);FD77:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X7);FD78:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X7);

(EI*sinh(mu*x)*c8*mu^3-cosh(mu*x)*L^3)/L^3

(EI*cosh(mu*x)*c8*mu^3-sinh(mu*x)*L^3)/L^3

(EI*sin(mu*x)*c8*mu^3-L^3*cos(mu*x))/L^3

(-EI*cos(mu*x)*c8*mu^3-sin(mu*x)*L^3)/L^3

(33)

>

FD81:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD82:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD83:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD84:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD85:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X8);FD86:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X8);FD87:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X8);FD88:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X8);

(-EI*cosh(mu*x)*c44*mu^2-L*sinh(mu*x)*mu)/L^2

(-EI*sinh(mu*x)*c44*mu^2-L*cosh(mu*x)*mu)/L^2

(EI*cos(mu*x)*c44*mu^2+L*sin(mu*x)*mu)/L^2

(EI*sin(mu*x)*c44*mu^2-L*cos(mu*x)*mu)/L^2

(34)

>

MM:=matrix(8,8,[[FD11,FD12,FD13,FD14,FD15,FD16,FD17,FD18],[FD21,FD22,FD23,FD24,FD25,FD26,FD27,FD28],[FD31,FD32,FD33,FD34,FD35,FD36,FD37,FD38],[FD41,FD42,FD43,FD44,FD45,FD46,FD47,FD48],[FD51,FD52,FD53,FD54,FD55,FD56,FD57,FD58],[FD61,FD62,FD63,FD64,FD65,FD66,FD67,FD68],[FD71,FD72,FD73,FD74,FD75,FD76,FD77,FD78],[FD81,FD82,FD83,FD84,FD85,FD86,FD87,FD88]]);

$MM := Matrix(8, 8, {(1, 1) = FD11, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 1, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = cosh(mu), (3, 6) = sinh(mu), (3, 7) = -cos(mu), (3, 8) = -sin(mu), (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = -sinh(mu), (4, 6) = -cosh(mu), (4, 7) = -sin(mu), (4, 8) = cos(mu), (5, 1) = cosh(mu*x), (5, 2) = sinh(mu*x), (5, 3) = -cos(mu*x), (5, 4) = -sin(mu*x), (5, 5) = -cosh(mu*x), (5, 6) = -sinh(mu*x), (5, 7) = cos(mu*x), (5, 8) = sin(mu*x), (6, 1) = -sinh(mu*x), (6, 2) = -cosh(mu*x), (6, 3) = -sin(mu*x), (6, 4) = cos(mu*x), (6, 5) = sinh(mu*x), (6, 6) = cosh(mu*x), (6, 7) = sin(mu*x), (6, 8) = -cos(mu*x), (7, 1) = (EI*sinh(mu*x)*c8*mu^3-cosh(mu*x)*L^3)/L^3, (7, 2) = (EI*cosh(mu*x)*c8*mu^3-sinh(mu*x)*L^3)/L^3, (7, 3) = (EI*sin(mu*x)*c8*mu^3-L^3*cos(mu*x))/L^3, (7, 4) = (-EI*cos(mu*x)*c8*mu^3-sin(mu*x)*L^3)/L^3, (7, 5) = cosh(mu*x), (7, 6) = sinh(mu*x), (7, 7) = cos(mu*x), (7, 8) = sin(mu*x), (8, 1) = (-EI*cosh(mu*x)*c44*mu^2-L*sinh(mu*x)*mu)/L^2, (8, 2) = (-EI*sinh(mu*x)*c44*mu^2-L*cosh(mu*x)*mu)/L^2, (8, 3) = (EI*cos(mu*x)*c44*mu^2+L*sin(mu*x)*mu)/L^2, (8, 4) = (EI*sin(mu*x)*c44*mu^2-L*cos(mu*x)*mu)/L^2, (8, 5) = sinh(mu*x)*mu/L, (8, 6) = cosh(mu*x)*mu/L, (8, 7) = -sin(mu*x)*mu/L, (8, 8) = cos(mu*x)*mu/L})$

(35)

Program end

>

Download Vibration_of_a_cracked_composite_beam.mw

>	restart: with(LinearAlgebra):

>	# Motion equation ( Vibration of a cracked composite beam using general solution) Edited by Adjal Yassine #

>	####################################################################

Motion equation of flexural vibration in normalized form

>	EIW^(iv)-momega^2*W=0;

(1)

The general solution form of bending vibration equation

>	W1:=A[1]cosh(mux)+A[2]sinh(mux)+A[3]cos(mux)+A[4]sin(mux);

(2)

where

>	E:=2682e6;L:=0.18;h:=0.004;b:=0.02;rho:=2600;area=bh;m:=rhohb;II:=(hb^3)/12:

(3)

>	mu:=((momega^2L^4/EI)^(1/4)):

Expression of cross-sectional rotation , the bending moment shear force and torsional moment are given as follows respectively

>	theta1 := (1/L)(A[1]musinh(mux)+A[2]mucosh(mux)-A[3]musin(mux)+A[4]mucos(mu*x));

(A[1]*mu*sinh(mu*x)+A[2]*mu*cosh(mu*x)-A[3]*mu*sin(mu*x)+A[4]*mu*cos(mu*x))/L

(4)

>	M1:= (EI/L^2)(A[1]mu^2cosh(mux)+A[2]mu^2sinh(mux)-A[3]mu^2cos(mux)-A[4]mu^2sin(mu*x));

EI*(A[1]*mu^2*cosh(mu*x)+A[2]*mu^2*sinh(mu*x)-A[3]*mu^2*cos(mu*x)-A[4]*mu^2*sin(mu*x))/L^2

(5)

>	S1:= (-EI/L^3)(A[1]mu^3sinh(mux)+A[2]mu^3cosh(mux)+A[3]mu^3sin(mux)-A[4]mu^3cos(mu*x));

-EI*(A[1]*mu^3*sinh(mu*x)+A[2]*mu^3*cosh(mu*x)+A[3]*mu^3*sin(mu*x)-A[4]*mu^3*cos(mu*x))/L^3

(6)

>

>	W2:=A[5]cosh(mux)+A[6]sinh(mux)+A[7]cos(mux)+A[8]sin(mux);

(7)

>

>	theta2:= (1/L)(A[5]musinh(mux)+A[6]mucosh(mux)-A[7]musin(mux)+A[8]mucos(mu*x));

(A[5]*mu*sinh(mu*x)+A[6]*mu*cosh(mu*x)-A[7]*mu*sin(mu*x)+A[8]*mu*cos(mu*x))/L

(8)

>	M2:= (EI/L^2)(A[5]mu^2cosh(mux)+A[6]mu^2sinh(mux)-A[7]mu^2cos(mux)-A[8]mu^2sin(mu*x));

EI*(A[5]*mu^2*cosh(mu*x)+A[6]*mu^2*sinh(mu*x)-A[7]*mu^2*cos(mu*x)-A[8]*mu^2*sin(mu*x))/L^2

(9)

>	S2:= -(EI/L^3)(A[5]mu^3sinh(mux)+A[6]mu^3cosh(mux)+A[7]mu^3sin(mux)-A[8]mu^3cos(mu*x));

-EI*(A[5]*mu^3*sinh(mu*x)+A[6]*mu^3*cosh(mu*x)+A[7]*mu^3*sin(mu*x)-A[8]*mu^3*cos(mu*x))/L^3

(10)

>

The boundary conditions at fixed end W1(0)=Theta(0)=0

>	X1:=eval(subs(x=0,W1));

(11)

>	X2:=eval(subs(x=0,theta1));

(mu*A[2]+mu*A[4])/L

(12)

>	X2:=collect(X2,mu)*(L/mu);

(13)

>

The boundary condtions at free end M2(1)=S2(1)=0

>	X3:=eval(subs(x=1,M2));

EI*(A[5]*mu^2*cosh(mu)+A[6]*mu^2*sinh(mu)-A[7]*mu^2*cos(mu)-A[8]*mu^2*sin(mu))/L^2

(14)

>	X3:=collect(X3,mu)*(L^2/mu^2/EI);

(15)

>	X4:=eval(subs(x=1,S2));

-EI*(A[5]*mu^3*sinh(mu)+A[6]*mu^3*cosh(mu)+A[7]*mu^3*sin(mu)-A[8]*mu^3*cos(mu))/L^3

(16)

>	X4:=collect(X4,mu);

-EI*(cosh(mu)*A[6]+sinh(mu)*A[5]-cos(mu)*A[8]+sin(mu)*A[7])*mu^3/L^3

(17)

>	X4:=collect(X4,EI)*(L^3/mu^3/EI);

(18)

>

The additional boundary conditions at crack location

>	X5:=combine(M1-M2);

(EI*cosh(mu*x)*mu^2*A[1]-EI*cosh(mu*x)*mu^2*A[5]+EI*sinh(mu*x)*mu^2*A[2]-EI*sinh(mu*x)*mu^2*A[6]-EI*cos(mu*x)*mu^2*A[3]+EI*cos(mu*x)*mu^2*A[7]-EI*sin(mu*x)*mu^2*A[4]+EI*sin(mu*x)*mu^2*A[8])/L^2

(19)

>	X5:=collect(X5,mu);

(EI*cosh(mu*x)*A[1]-EI*cosh(mu*x)*A[5]+EI*sinh(mu*x)*A[2]-EI*sinh(mu*x)*A[6]-cos(mu*x)*EI*A[3]+A[7]*cos(mu*x)*EI-A[4]*sin(mu*x)*EI+A[8]*sin(mu*x)*EI)*mu^2/L^2

(20)

>	X5:=collect(X5,EI)*(L^2/mu^2/EI);

A[1]*cosh(mu*x)-A[5]*cosh(mu*x)+A[2]*sinh(mu*x)-A[6]*sinh(mu*x)-A[3]*cos(mu*x)+A[7]*cos(mu*x)-A[4]*sin(mu*x)+A[8]*sin(mu*x)

(21)

>	X6:=combine(S1-S2);

(-EI*cosh(mu*x)*mu^3*A[2]+EI*cosh(mu*x)*mu^3*A[6]-EI*sinh(mu*x)*mu^3*A[1]+EI*sinh(mu*x)*mu^3*A[5]+EI*cos(mu*x)*mu^3*A[4]-EI*cos(mu*x)*mu^3*A[8]-EI*sin(mu*x)*mu^3*A[3]+EI*sin(mu*x)*mu^3*A[7])/L^3

(22)

>	X6:=collect(X6,mu);

(-EI*cosh(mu*x)*A[2]+EI*cosh(mu*x)*A[6]-EI*sinh(mu*x)*A[1]+EI*A[5]*sinh(mu*x)+cos(mu*x)*A[4]*EI-cos(mu*x)*A[8]*EI-sin(mu*x)*EI*A[3]+sin(mu*x)*A[7]*EI)*mu^3/L^3

(23)

>	X6:=collect(X6,EI)*(L^3/mu^3/EI);

-cosh(mu*x)*A[2]+cosh(mu*x)*A[6]-sinh(mu*x)*A[1]+sinh(mu*x)*A[5]+cos(mu*x)*A[4]-cos(mu*x)*A[8]-sin(mu*x)*A[3]+sin(mu*x)*A[7]

(24)

>

>	X7:=combine(W2-(W1+c8*S1));

(EI*cosh(mu*x)*c8*mu^3*A[2]+EI*sinh(mu*x)*c8*mu^3*A[1]-EI*cos(mu*x)*c8*mu^3*A[4]+EI*sin(mu*x)*c8*mu^3*A[3]-cosh(mu*x)*A[1]*L^3+cosh(mu*x)*A[5]*L^3-sinh(mu*x)*A[2]*L^3+sinh(mu*x)*A[6]*L^3-cos(mu*x)*A[3]*L^3+cos(mu*x)*A[7]*L^3-sin(mu*x)*A[4]*L^3+sin(mu*x)*A[8]*L^3)/L^3

(25)

>	X8:=combine (theta2-(theta1+c44*M1));

(-EI*cosh(mu*x)*c44*mu^2*A[1]-EI*sinh(mu*x)*c44*mu^2*A[2]+EI*cos(mu*x)*c44*mu^2*A[3]+EI*sin(mu*x)*c44*mu^2*A[4]-L*cosh(mu*x)*mu*A[2]+L*cosh(mu*x)*mu*A[6]-L*sinh(mu*x)*mu*A[1]+L*sinh(mu*x)*mu*A[5]-L*cos(mu*x)*mu*A[4]+L*cos(mu*x)*mu*A[8]+L*sin(mu*x)*mu*A[3]-L*sin(mu*x)*mu*A[7])/L^2

(26)

>

The characteristic matrix function of frequency

>	FD8:=subs(A[1]=1,A[3]=0,X1);FD12:=subs(A[1]=0,A[3]=0,X1);FD13:=subs(A[1]=0,A[3]=1,X1);FD14:=subs(A[1]=0,A[3]=0,X1);FD15:=subs(A[1]=0,A[3]=0,X1);FD16:=subs(A[1]=0,A[3]=0,X1);FD17:=subs(A[1]=0,A[3]=0,X1);FD18:=subs(A[1]=0,A[3]=0,X1);

(27)

>	FD21:=subs(A[2]=0,A[4]=0,X2);FD22:=subs(A[2]=1,A[4]=0,X2);FD23:=subs(A[2]=0,A[4]=0,X2);FD24:=subs(A[2]=0,A[4]=1,X2);FD25:=subs(A[2]=0,A[4]=0,X2);FD26:=subs(A[2]=0,A[4]=0,X2);FD27:=subs(A[2]=0,A[4]=0,X2);FD28:=subs(A[2]=0,A[4]=0,X2);

(28)

>

FD31:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD32:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD33:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD34:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X3);FD35:=subs(A[5]=1,A[6]=0,A[7]=0,A[8]=0,X3);;FD36:=subs(A[5]=0,A[6]=1,A[7]=0,A[8]=0,X3);FD37:=subs(A[5]=0,A[6]=0,A[7]=1,A[8]=0,X3);FD38:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=1,X3);

(29)

>

FD41:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD42:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD43:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD44:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=0,X4);FD45:=subs(A[5]=1,A[6]=0,A[7]=0,A[8]=0,X4);FD46:=subs(A[5]=0,A[6]=1,A[7]=0,A[8]=0,X4);FD47:=subs(A[5]=0,A[6]=0,A[7]=1,A[8]=0,X4);FD48:=subs(A[5]=0,A[6]=0,A[7]=0,A[8]=1,X4);

(30)

>

FD51:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD52:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD53:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD54:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X5);FD55:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X5);FD56:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X5);FD57:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X5);FD58:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X5);

(31)

>

FD61:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD62:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD63:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD64:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X6);FD65:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X6);FD66:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X6);FD67:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X6);FD68:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X6);

(32)

>

FD71:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD72:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD73:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD74:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X7);FD75:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X7);FD76:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X7);FD77:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X7);FD78:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X7);

(EI*sinh(mu*x)*c8*mu^3-cosh(mu*x)*L^3)/L^3

(EI*cosh(mu*x)*c8*mu^3-sinh(mu*x)*L^3)/L^3

(EI*sin(mu*x)*c8*mu^3-L^3*cos(mu*x))/L^3

(-EI*cos(mu*x)*c8*mu^3-sin(mu*x)*L^3)/L^3

(33)

>

FD81:=subs(A[1]=1,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD82:=subs(A[1]=0,A[2]=1,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD83:=subs(A[1]=0,A[2]=0,A[3]=1,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD84:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=1,A[5]=0,A[6]=0,A[7]=0,A[8]=0,X8);FD85:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=1,A[6]=0,A[7]=0,A[8]=0,X8);FD86:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=1,A[7]=0,A[8]=0,X8);FD87:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=1,A[8]=0,X8);FD88:=subs(A[1]=0,A[2]=0,A[3]=0,A[4]=0,A[5]=0,A[6]=0,A[7]=0,A[8]=1,X8);

(-EI*cosh(mu*x)*c44*mu^2-L*sinh(mu*x)*mu)/L^2

(-EI*sinh(mu*x)*c44*mu^2-L*cosh(mu*x)*mu)/L^2

(EI*cos(mu*x)*c44*mu^2+L*sin(mu*x)*mu)/L^2

(EI*sin(mu*x)*c44*mu^2-L*cos(mu*x)*mu)/L^2

(34)

>

MM:=matrix(8,8,[[FD11,FD12,FD13,FD14,FD15,FD16,FD17,FD18],[FD21,FD22,FD23,FD24,FD25,FD26,FD27,FD28],[FD31,FD32,FD33,FD34,FD35,FD36,FD37,FD38],[FD41,FD42,FD43,FD44,FD45,FD46,FD47,FD48],[FD51,FD52,FD53,FD54,FD55,FD56,FD57,FD58],[FD61,FD62,FD63,FD64,FD65,FD66,FD67,FD68],[FD71,FD72,FD73,FD74,FD75,FD76,FD77,FD78],[FD81,FD82,FD83,FD84,FD85,FD86,FD87,FD88]]);

$MM := Matrix(8, 8, {(1, 1) = FD11, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 1, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = cosh(mu), (3, 6) = sinh(mu), (3, 7) = -cos(mu), (3, 8) = -sin(mu), (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = -sinh(mu), (4, 6) = -cosh(mu), (4, 7) = -sin(mu), (4, 8) = cos(mu), (5, 1) = cosh(mu*x), (5, 2) = sinh(mu*x), (5, 3) = -cos(mu*x), (5, 4) = -sin(mu*x), (5, 5) = -cosh(mu*x), (5, 6) = -sinh(mu*x), (5, 7) = cos(mu*x), (5, 8) = sin(mu*x), (6, 1) = -sinh(mu*x), (6, 2) = -cosh(mu*x), (6, 3) = -sin(mu*x), (6, 4) = cos(mu*x), (6, 5) = sinh(mu*x), (6, 6) = cosh(mu*x), (6, 7) = sin(mu*x), (6, 8) = -cos(mu*x), (7, 1) = (EI*sinh(mu*x)*c8*mu^3-cosh(mu*x)*L^3)/L^3, (7, 2) = (EI*cosh(mu*x)*c8*mu^3-sinh(mu*x)*L^3)/L^3, (7, 3) = (EI*sin(mu*x)*c8*mu^3-L^3*cos(mu*x))/L^3, (7, 4) = (-EI*cos(mu*x)*c8*mu^3-sin(mu*x)*L^3)/L^3, (7, 5) = cosh(mu*x), (7, 6) = sinh(mu*x), (7, 7) = cos(mu*x), (7, 8) = sin(mu*x), (8, 1) = (-EI*cosh(mu*x)*c44*mu^2-L*sinh(mu*x)*mu)/L^2, (8, 2) = (-EI*sinh(mu*x)*c44*mu^2-L*cosh(mu*x)*mu)/L^2, (8, 3) = (EI*cos(mu*x)*c44*mu^2+L*sin(mu*x)*mu)/L^2, (8, 4) = (EI*sin(mu*x)*c44*mu^2-L*cos(mu*x)*mu)/L^2, (8, 5) = sinh(mu*x)*mu/L, (8, 6) = cosh(mu*x)*mu/L, (8, 7) = -sin(mu*x)*mu/L, (8, 8) = cos(mu*x)*mu/L})$

(35)

Program end

>

Download Vibration_of_a_cracked_composite_beam.mwVibration_of_a_cracked_composite_beam.mwVibration_of_a_cracked_composite_beam.mw

Splitting PDE parameterized symmetries and parameter...

by: ecterrab 14540 Product: Maple

October 18 2019

5 6

Splitting PDE parameterized symmetries

and Parameter-continuous symmetry transformations

The determination of symmetries for partial differential equation systems (PDE) is relevant in several contexts, the most obvious of which is of course the determination of the PDE solutions. For instance, generally speaking, the knowledge of a N-dimensional Lie symmetry group can be used to reduce the number of independent variables of PDE by N. So if PDE depends only on N independent variables, that amounts to completely solving it. If only N-1 symmetries are known or can be successfully used then PDE becomes and ODE; etc., all advantageous situations. In Maple, a complete set of symmetry commands, to perform each step of the symmetry approach or several of them in one go, is part of the PDEtools package.

Besides the dependent and independent variables, PDE frequently depends on some constant parameters, and besides the PDE symmetries for arbitrary values of those parameters, for some particular values of them, PDE transforms into a completely different problem, admitting different symmetries. The question then is: how can you determine those particular values of the parameters and the corresponding different symmetries? That was the underlying subject of a recent question in Mapleprimes. The answer to those questions is relatively simple and yet not entirely obvious for most of us, motivating this post, organized briefly around one example.

To reproduce the input/output below you need Maple 2019 and to have installed the Physics Updates v.449 or higher.

Consider the family of Korteweg-de Vries equation for involving three constant parameters . For convenience (simpler input and more readable output) use the diff_table and declare commands

>

(1)

>

(2)

This pde admits a 4-dimensional symmetry group, whose infinitesimals - for arbitrary values of the parameters - are given by

>

(3)

Looking at pde (1) as a nonlinear problem in u, a, b and q, it splits into four cases for some particular values of the parameter:

>	$pde__cases := casesplit(bu(x, t)(diff(u(x, t), x))+a(diff(u(x, t), x))+q(diff(diff(diff(u(x, t), x), x), x))+diff(u(x, t), t) = 0, parameters = {a, b, q}, caseplot)$

`casesplit/ans`([diff(diff(diff(u(x, t), x), x), x) = -(b*u(x, t)*(diff(u(x, t), x))+a*(diff(u(x, t), x))+diff(u(x, t), t))/q], [q <> 0]), `casesplit/ans`([diff(u(x, t), x) = -(diff(u(x, t), t))/(b*u(x, t)+a), q = 0], [b*u(x, t)+a <> 0]), `casesplit/ans`([u(x, t) = -a/b, q = 0], [b <> 0]), `casesplit/ans`([diff(u(x, t), t) = 0, a = 0, b = 0, q = 0], [])

(4)

The legend above indicates the pivots and the tree of cases, depending on whether each pivot is equal or different from 0. At the end there is the algebraic sequence of cases. The first case is the general case, for which the symmetry infinitesimals were computed as above, but clearly the other three cases admit more general symmetries. Consider for instance the second case, pass the to ignore the parameterizing equation , and you get

>

[_xi[x](x, t, u) = _F3(x, t, u), _xi[t](x, t, u) = Intat(((b*u+a)*(D[1](_F3))(_a, ((b*u+a)*t-x+_a)/(b*u+a), u)-_F1(u, ((b*u+a)*t-x)/(b*u+a))*b+(D[2](_F3))(_a, ((b*u+a)*t-x+_a)/(b*u+a), u))/(b*u+a)^2, _a = x)+_F2(u, ((b*u+a)*t-x)/(b*u+a)), _eta[u](x, t, u) = _F1(u, ((b*u+a)*t-x)/(b*u+a))]

(5)

These infinitesimals are indeed much more general than , in fact so general that (5) is almost unreadable ... Specialize the three arbitrary functions into something "easy" just to be able follow - e.g. take _F1 to be just the + operator, _F2 the * operator and

>

[_xi[x](x, t, u) = 1, _xi[t](x, t, u) = Intat(-(u+((b*u+a)*t-x)/(b*u+a))*b/(b*u+a)^2, _a = x)+u*((b*u+a)*t-x)/(b*u+a), _eta[u](x, t, u) = u+((b*u+a)*t-x)/(b*u+a)]

(6)

>

[_xi[x](x, t, u) = 1, _xi[t](x, t, u) = (b^3*t*u^4+((3*a*t-x)*u^3-u^2*x-t*x*u)*b^2+((3*a^2*t-2*a*x)*u^2-a*u*x-a*t*x+x^2)*b+a^2*u*(a*t-x))/(b*u+a)^3, _eta[u](x, t, u) = (b*u^2+(b*t+a)*u+a*t-x)/(b*u+a)]

(7)

This symmetry is of course completely different than computed for the general case.

The symmetry (7) can be verified against or directly against pde after substituting .

(8)

>

(9)

>

(10)

Summarizing: "to split PDE symmetries into cases according to the values of the PDE parameters, split the PDE into cases with respect to these parameters (command PDEtools:-casesplit ) then compute the symmetries for each case"

Parameter continuous symmetry transformations

A different, however closely related question, is whether pde admits "symmetries with respect to the parameters a, b and q", so whether exists continuous transformations of the parameters a, b and q that leave pde invariant in form.

Beforehand, note that since the parameters are constants with regards to the dependent and independent variables (here ), such continuous symmetry transformations cannot be used directly to compute a solution for pde. They can, however, be used to reduce the number of parameters. And in some contexts, that is exactly what we need, for example to entirely remove the splitting into cases due to their presence, or to proceed applying a solving method that is valid only when there are no parameters (frequently the case when computing exact solutions to "PDE & Boundary Conditions").

To compute such "continuous symmetry transformations of the parameters" that leave pde invariant one can always think of these parameters as "additional independent variables of pde". In terms of formulation, that amounts to replacing the dependency in the dependent variable, i.e. replace by

>

(11)

Compute now the infinitesimals: note there are now three additional ones, related to continuous transformations of and - for readability, avoid displaying the redundant functionality on the left-hand sides of these infinitesimals

>

[_xi[x] = (1/3)*(_F4(a, b, q)*q+_F3(a, b, q))*x/q+_F6(a, b, q)*t+_F7(a, b, q), _xi[t] = _F4(a, b, q)*t+_F5(a, b, q), _xi[a] = _F1(a, b, q), _xi[b] = _F2(a, b, q), _xi[q] = _F3(a, b, q), _eta[u] = (1/3)*((b*u+a)*_F3(a, b, q)-2*((b*u+a)*_F4(a, b, q)+(3/2)*u*_F2(a, b, q)+(3/2)*_F1(a, b, q)-(3/2)*_F6(a, b, q))*q)/(b*q)]

(12)

This result is more general than what is convenient for algebraic manipulations, so specialize the seven arbitrary functions of and keep only the first symmetry that result from this specialization: that suffices to illustrate the removal of any of the three parameters a, b, or q

>

(13)

To remove the parameters, as it is standard in the symmetry approach, compute a transformation to canonical coordinates, with respect to the parameter a. That means a transformation that changes the list of infinitesimals, or likewise its infinitesimal generator representation,

>

proc (f) options operator, arrow; diff(f, a)-(diff(f, u))/b end proc

(14)

into or its equivalent generator representation

That same transformation, when applied to , entirely removes the parameter a.

The transformation is computed using CanonicalCoordinates and the last argument indicates the "independent variable" (in our case a parameter) that the transformation should remove. We choose to remove a

>

(15)

>

(16)

Invert this transformation in order to apply it

>

(17)

The next step is not necessary, but just to understand how all this works, verify its action over the infinitesimal generator proc (f) options operator, arrow; diff(f, a)-(diff(f, u))/b end proc

>

ChangeSymmetry({a = alpha, b = beta, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*beta-alpha)/beta}, proc (f) options operator, arrow; diff(f, a)-(diff(f, u))/b end proc, [upsilon(xi, tau, alpha, beta, chi), xi, tau, alpha, beta, chi])

(18)

Now that we see the transformation (17) is the one we want, just use it to change variables in

>	$PDEtools:-dchange({a = alpha, b = beta, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*beta-alpha)/beta}, pde__xtabq, [upsilon(xi, tau, alpha, beta, chi), xi, tau, alpha, beta, chi], simplify)$

(19)

As expected, this result depends only on two parameters, , and , and the one equivalent to a (that is , see the transformation used (17)), is not present anymore.

To remove b or q we use the same steps, (15), (17) and (19), just changing the parameter to be removed, indicated as the last argument in the call to CanonicalCoordinates . For example, to eliminate b (represented in the new variables by ), input

>

(20)

>

(21)

>	$PDEtools:-dchange({a = beta, b = alpha, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*alpha-beta)/alpha}, pde__xtabq, [upsilon(xi, tau, alpha, beta, chi), xi, tau, alpha, beta, chi], simplify)$

(22)

and as expected this result does not contain To remove a second parameter, the whole cycle is repeated starting with computing infinitesimals, for instance for (22). Finally, the case of function parameters is treated analogously, by considering the function parameters as additional dependent variables instead of independent ones.

Download How_to_split_symmetries_into_cases_(II).mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

Pole simulation with tensions in MapleSim

by: Lenin Araujo Castillo 1790 Products: Maple 2018 , MapleSim 2018

October 14 2019

2 0

Application developed using Maple and MapleSim. You can observe the vector analysis using Maple and the simulation using MapleSim. Also included a video of the result. It is a simple structure. A pole fastened by two cables and a force applied to the top. The results are to calculate tensions one and two. Consider the mass of each rope. In spanish.

POSTE_PARADO.zip

Lenin Araujo Castillo

Ambassador of Maple

Maple May benefit from starting a Stack Exchange...

by: Adam Ledger 360 Product: Maple

October 12 2019

6 7

I have noticed that there exists a Stack Exchange site for mathematica, and not for maple. My discussions with the part of Stack Exchange that handle the creation of a new Stack Exchange community have said that I must accrue a certain level of interest in the subject in order for it to be approved, and so I thought I would begin here to see if there is suffice level of interest.

This would not diminish the use of the Maple Primes forum, and an additional proposal, in consideration of the years of dedication that have gone into this domain, be to pool the data between the two, make reputation points the same on both, perhaps even user profiles and questions answered already linkable, and all of the questions already addressed here showing up in the search on both domains.

I am proposing this simply because I want to encourage the use of maple, and have noted that Stack Exchange is very popular.

So I am posting this to get overall feedback from other Maple users, as to what their opinion is regarding this proposal, and advice on whether it should and how it ought to be pursued.

Integral Transforms (revamped) and PDE

by: ecterrab 14540 Product: Maple

October 01 2019

8 14

Integral Transforms (revamped) and PDEs

Integral transforms, implemented in Maple as the inttrans package, are special integrals that appear frequently in mathematical-physics and that have remarkable properties. One of the main uses of integral transforms is for the computation of exact solutions to ordinary and partial differential equations with initial/boundary conditions. In Maple, that functionality is implemented in dsolve/inttrans and in pdsolve/boundary conditions .

During the last months, we have been working heavily on several aspects of these integral transform functions and this post is about that. This is work in progress, in collaboration with Katherina von Bulow.

The integral transforms are represented by the commands of the inttrans package:

>

(1)

Three of these commands, addtable, savetable, and setup (this one is new, only present after installing the Physics Updates) are "administrative" commands while the others are computational representations for integrals. For example,

>

[fourier(a, b, z) = Int(a/exp(I*b*z), b = -infinity .. infinity), MathematicalFunctions:-`with no restrictions on `(a, b, z)]

(2)

>

[mellin(a, b, z) = Int(a*b^(z-1), b = 0 .. infinity), MathematicalFunctions:-`with no restrictions on `(a, b, z)]

(3)

For all the integral transform commands, the first argument is the integrand, the second one is the dummy integration variable of a definite integral and the third one is the evaluation point. (also called transform variable). The integral representation is also visible using the convert network

>

laplace(f(t), t, s) = Int(f(t)*exp(-s*t), t = 0 .. infinity)

(4)

Having in mind the applications of these integral transforms to compute integrals and exact solutions to PDE with boundary conditions, five different aspects of these transforms received further development:

•	Compute Derivatives: Yes or No

•	Numerical Evaluation

•	Two Hankel Transform Definitions

•	More integral transform results

•	Mellin and Hankel transform solutions for Partial Differential Equations with boundary conditions

The project includes having all these tranforms available at user level (not ready), say as FourierTransform for inttrans:-fourier, so that we don't need to input anymore. Related to these changes we also intend to have not return anymore, and return itself instead, unevaluated, so that one can set its value according to the problem/preferred convention (typically 0, 1/2 or 1) and have all the Maple library following that choice.

The material presented in the following sections is reproducible already in Maple 2019 by installing the latest Physics Updates (v.435 or higher),

Compute derivatives: Yes or No.

For historical reasons, previous implementations of these integral transform commands did not follow a standard paradigm of computer algebra: "Given a function , the input should return the derivative of ". The implementation instead worked in the opposite direction: if you were to input the result of the derivative, you would receive the derivative representation. For example, to the input you would receive . This is particularly useful for the purpose of using integral transforms to solve differential equations but it is counter-intuitive and misleading; Maple knows the differentiation rule of these functions, but that rule was not evident anywhere. It was not clear how to compute the derivative (unless you knew the result in advance).

To solve this issue, a new command, setup, has been added to the package, so that you can set "whether or not" to compute derivatives, and the default has been changed to computederivatives = true while the old behavior is obtained only if you input . For example, after having installed the Physics Updates,

>

`The "Physics Updates" version in the MapleCloud is 435 and is the same as the version installed in this computer, created 2019, October 1, 12:46 hours, found in the directory /Users/ecterrab/maple/toolbox/2019/Physics Updates/lib/`

(1.1)

the current settings can be queried via

>

(1.2)

and so differentiating returns the derivative computed

>

(1.3)

while changing this setting to work as in previous releases you have this computation reversed: you input the output (1.3) and you get the corresponding input

>

(1.4)

>

(1.5)

Reset the value of computederivatives

>

(1.6)

>

(1.7)

In summary: by default, derivatives of integral transforms are now computed; if you need to work with these derivatives as in previous releases, you can input . This setting can be changed any time you want within one and the same Maple session, and changing it does not have any impact on the performance of intsolve, dsolve and pdsolve to solve differential equations using integral transforms.

>

Numerical Evaluation

In previous releases, integral transforms had no numerical evaluation implemented. This is in the process of changing. So, for example, to numerically evaluate the inverse laplace transform ( invlaplace command), three different algorithms have been implemented: Gaver-Stehfest, Talbot and Euler, following the presentation by Abate and Whitt, "Unified Framework for Numerically Inverting Laplace Transforms", INFORMS Journal on Computing 18(4), pp. 408–421, 2006.

For example, consider the exact solution to this partial differential equation subject to initial and boundary conditions

>

Note that these two conditions are not entirely compatible: the solution returned cannot be valid for and simultaneously. However, a solution discarding that point does exist and is given by

>

u(x, t) = -invlaplace(exp(-(1/2)*s^(1/2)*t)/s, s, x)+1

(2.1)

Verifying the solution, one condition remains to be tested

>

[0, 0, -invlaplace(exp(-(1/2)*s^(1/2)*t)/s, s, 0)]

(2.2)

Since we now have numerical evaluation rules, we can test that what looks different from 0 in the above is actually 0.

>

-invlaplace(exp(-(1/2)*s^(1/2)*t)/s, s, 0)

(2.3)

Add a small number to the initial value of t to skip the point

>

The default method used is the method of Euler sums and the numerical evaluation is performed as usual using the evalf command. For example, consider

>

The Laplace transform of F is given by

>

(1/2)*2^(1/2)*Pi^(1/2)*exp(-(1/2)/s)/s^(3/2)

(2.4)

and the inverse Laplace transform of LT in inert form is

>

%invlaplace((1/2)*2^(1/2)*Pi^(1/2)*exp(-(1/2)/s)/s^(3/2), s, t)

(2.5)

At we have

>

%invlaplace((1/2)*2^(1/2)*Pi^(1/2)*exp(-(1/2)/s)/s^(3/2), s, 1)

(2.6)

>

(2.7)

This result is consistent with the one we get if we first compute the exact form of the inverse Laplace transform at t = 1:

>

%invlaplace((1/2)*2^(1/2)*Pi^(1/2)*exp(-(1/2)/s)/s^(3/2), s, 1) = sin(2^(1/2))

(2.8)

>

(2.9)

In addition to the standard use of evalf to numerically evaluate inverse Laplace transforms, one can invoke each of the three different methods implemented using the MathematicalFunctions:-Evalf command

>

(2.10)

>

(2.11)

>

(2.12)

>

(2.13)

Regarding the method we use by default: from a numerical experiment with varied problems we have concluded that our implementation of the Euler (sums) method is faster and more accurate than the other two.

Two Hankel transform definitions

In previous Maple releases, the definition of the Hankel transform was given by

hankel(f(t), t, s, nu) = Int(f(t)*sqrt(s*t)*BesselJ(nu, s*t), t = 0 .. infinity)

where is the function. This definition, sometimes called alternative definition of the Hankel transform, has the inconvenience of the square root in the integrand, complicating the form of the hankel transform for the Laplacian in cylindrical coordinates. On the other hand, the definition more frequently used in the literature is

hankel(f(t), t, s, nu) = Int(f(t)*t*BesselJ(nu, s*t), t = 0 .. infinity)

With it, the Hankel transform of diff(u(r, t), r, r)+(diff(u(r, t), r))/r+diff(u(r, t), t, t) is given by the simple ODE form d^2*`&Hopf;`(k, t)/dt^2-k^2*`&Hopf;`(k, t) . Not just this transform but several other ones acquire a simpler form with the definition that does not have a square root in the integrand.

So the idea is to align Maple with this simpler definition, while keeping the previous definition as an alternative. Hence, by default, when you load the inttrans package, the new definition in use for the Hankel transform is

>

hankel(f(t), t, s, nu) = Int(f(t)*t*BesselJ(nu, s*t), t = 0 .. infinity)

(3.1)

You can change this default so that Maple works with the alternative definition as in previous releases. For that purpose, use the new inttrans:-setup command (which you can also use to query about the definition in use at any moment):

>

(3.2)

This change in definition is automatically taken into account by other parts of the Maple library using the Hankel transform. For example, the differentiation rule with the new definition is

>

(3.3)

This differentiation rule resembles (is connected to) the differentiation rule for BesselJ, and this is another advantage of the new definition.

>

%diff(BesselJ(nu, z), z) = -BesselJ(nu+1, z)+nu*BesselJ(nu, z)/z

(3.4)

Furthermore, several transforms have acquired a simpler form, as for example:

>

%hankel(exp(I*a*r)/r, r, k, 0) = 1/(-a^2+k^2)^(1/2)

(3.5)

Let's compare: make the definition be as in previous releases.

>

(3.6)

>

hankel(f(t), t, s, nu) = Int(f(t)*(s*t)^(1/2)*BesselJ(nu, s*t), t = 0 .. infinity)

(3.7)

The differentiation rule with the previous (alternative) definition was not as simple:

>

(3.8)

And the transform (3.5) was also not so simple:

>

%hankel(exp(I*a*r)/r, r, k, 0) = (I*a*hypergeom([3/4, 3/4], [3/2], a^2/k^2)*GAMMA(3/4)^4+Pi^2*k*hypergeom([1/4, 1/4], [1/2], a^2/k^2))/(k*Pi*GAMMA(3/4)^2)

(3.9)

Reset to the new default value of the definition.

>

(3.10)

>

hankel(f(t), t, s, nu) = Int(f(t)*t*BesselJ(nu, s*t), t = 0 .. infinity)

(3.11)

More integral transform results

The revision of the integral transforms includes also filling gaps: previous transforms that were not being computed are now computed. Still with the Hankel transform, consider the operators

>

`D/t` := proc (u) options operator, arrow; (diff(u, t))/t end proc

Being able to transform these operators into algebraic expressions or differential equations of lower order is key for solving PDE problems with Boundary Conditions.

Consider, for instance, this ODE

>

(4.1)

>

((diff(diff(diff(u(t), t), t), t))*t^3-12*(diff(diff(u(t), t), t))*t^2+57*(diff(u(t), t))*t-105*u(t))/t^3

(4.2)

Its Hankel transform is a simple algebraic expression

>

(4.3)

An example with formula_plus

>

((diff(diff(diff(diff(u(t), t), t), t), t))*t^4+38*(diff(diff(diff(u(t), t), t), t))*t^3+477*(diff(diff(u(t), t), t))*t^2+2295*(diff(u(t), t))*t+3465*u(t))/t^4

(4.4)

>

(4.5)

In the case of hankel , not just differential operators but also several new transforms are now computable

>

piecewise(nu = 0, Dirac(k)/k, nu/k^2)

(4.6)

>

piecewise(And(nu = 0, m = 0), Dirac(k)/k, 2^(m+1)*k^(-m-2)*GAMMA(1+(1/2)*m+(1/2)*nu)/GAMMA((1/2)*nu-(1/2)*m))

(4.7)

>

Mellin and Hankel transform solutions for Partial Differential Equations with Boundary Conditions

In previous Maple releases, the Fourier and Laplace transforms were used to compute exact solutions to PDE problems with boundary conditions. Now, Mellin and Hankel transforms are also used for that same purpose.

Example:

>

pde := x^2*(diff(u(x, y), x, x))+x*(diff(u(x, y), x))+diff(u(x, y), y, y) = 0

iv := u(x, 0) = 0, u(x, 1) = piecewise(0 <= x and x < 1, 1, 1 < x, 0)

>

(5.1)

As usual, you can let pdsolve choose the solving method, or indicate the method yourself:

>

pde := diff(u(r, t), r, r)+(diff(u(r, t), r))/r+diff(u(r, t), t, t) = -Q__0*q(r)
iv := u(r, 0) = 0

>

u(r, t) = Q__0*(-hankel(exp(-s*t)*hankel(q(r), r, s, 0)/s^2, s, r, 0)+hankel(hankel(q(r), r, s, 0)/s^2, s, r, 0))

(5.2)

It is sometimes preferable to see these solutions in terms of more familiar integrals. For that purpose, use

>

u(r, t) = Q__0*(-(Int(exp(-s*t)*(Int(q(r)*r*BesselJ(0, r*s), r = 0 .. infinity))*BesselJ(0, r*s)/s, s = 0 .. infinity))+Int((Int(q(r)*r*BesselJ(0, r*s), r = 0 .. infinity))*BesselJ(0, r*s)/s, s = 0 .. infinity))

(5.3)

An example where the hankel transform is computable to the end:

>

pde := c^2*(diff(u(r, t), r, r)+(diff(u(r, t), r))/r) = diff(u(r, t), t, t)
iv := u(r, 0) = A*a/(a^2+r^2)^(1/2), (D[2](u))(r, 0) = 0

>

u(r, t) = (1/2)*A*a*((-c^2*t^2+(2*I)*a*c*t+a^2+r^2)^(1/2)+(-c^2*t^2-(2*I)*a*c*t+a^2+r^2)^(1/2))/((-c^2*t^2-(2*I)*a*c*t+a^2+r^2)^(1/2)*(-c^2*t^2+(2*I)*a*c*t+a^2+r^2)^(1/2))

(5.4)

>

Download Integral_Transforms_(revamped).mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

Graph Theory and Pokémon Application

by: vdbustos 40 Product: Maple 2019

September 19 2019

8 0

Hi there!

One of my favorite videogames is pokémon as you can probably guess from the title. As a player I always wanted to optimize my chances of obtaining the rarest and best pokémon in the game. I have been working on an application that aims to use graph theory to analyze the game Pokémon Blue. The application explores the following questions:

Which is the rarest pokémon in the game?
Where can I find an specific pokémon and with what probabilities?
What is the place with most different species of wild pokémon?

I also included algorithms for the following: Given a certain desired team

Find the minimum amount of places to visit to catch them and return the list of the places the player will need to visit.
What are the routes with best probabilities to catch each pokémon from my desired team?

Check out my application at: https://www.maplesoft.com/applications/view.aspx?SID=154565.

The following are some of the results obtained in the app:

What is the most common pokémon?

I did not only considered the amount of places a pokémon can appear in but also the probabilities of it appearing in each place.

What are the connections between pokémon and places?

In my graph, I connected a pokémon and a place if such pokémon could be caught in that place. The following is an example for the pokémon Pidgey. The weights of the edges are the probabilities of finding Pidgey in each route.

Viceversa, I did the same for how a route is connected to the pokémons in it:

Map of the Game
I also generated a colour coded version for the map of the game: where blue means that the place is a water route, brown means it's a cave and green means it's a tall-grass route.
It's amazing what Maple's graph theory toolbox can do.

Arabic Cipher

by: Fereydoon_Shekofte 223 Product: Maple 18

September 15 2019

3 11

This is a simple encryption method to hide text messages

Mentioned in Arabic manuscrips with more than hundreds years old ...

PRINCIPLE :

Just the place of letters in the sentence rearranged as described below :

For example "ABCDE" we pick up the First letter "A" from the left and write it as the last letter in the Right "......A"

but this time we pick up the letter "E" as the last letter from Right and place it at the Left Side of the previous one ".....EA"

and this cycle continue until for rest letters ... "CDBEA" .

by this way the text become hard to discover !

It is Amazing that for decoding this message you should repeat the same rearrangment algorithm several times until the readable text appears as the first "ABCDE"

EXample :

"AlbertEinstein"

"iEntsrteebilnA"

"eterbsitlnnEAi"

"tilsnbnrEeAtie"

"rnEbenAstliiet"

"sAtnleibiEentr"

"biieElennttArs"

"nenltEteAirisb"

"etAEitrlinsebn"

"lritnisEeAbtne"

"EseiAnbttinrel"

"tbtniAnireeslE"

"inrAeienstlbEt"

"nesitelAbrEnti"

"AlbertEinstein"

the same text appeared after 14 step cycle

Arabic Cipher

>

ArabicCipher := proc (x) options operator, arrow; StringTools[Permute](x, [seq(1+iquo(StringTools[Length](x), 2)+((1/2)*i+(1/2)*irem(i, 2))*(-1)^(i+irem(StringTools[Length](x), 2)), i = 0 .. StringTools[Length](x)-1)]) end proc

proc (x) options operator, arrow; StringTools[Permute](x, [seq(1+iquo(StringTools[Length](x), 2)+((1/2)*i+(1/2)*irem(i, 2))*(-1)^(i+irem(StringTools[Length](x), 2)), i = 0 .. StringTools[Length](x)-1)]) end proc

(1.1)

>

"iEntsrteebilnA", "eterbsitlnnEAi", "tilsnbnrEeAtie", "rnEbenAstliiet", "sAtnleibiEentr", "biieElennttArs", "nenltEteAirisb", "etAEitrlinsebn", "lritnisEeAbtne", "EseiAnbttinrel", "tbtniAnireeslE", "inrAeienstlbEt", "nesitelAbrEnti", "AlbertEinstein"

(1.2)

>

"nSohoedkyoefrteeF", "yokedferotheoeSFn", "otrheefodeeSkFony", "deoefSekeFhorntyo", "eFkheoSrfnetoyeod", "fnreStooeyhekoFde", "eyohoetkSoeFrdnef", "SoketFerodhnoeyfe", "odrhenFoteeykfoeS", "teoeFynkefhoredSo", "efkhnoyrFeedoSeot", "FereydoonShekofte"

(1.3)

>

Download Arabic_Cipher.mw

Introducing the Maple Companion app!

by: Karishma 805 Product: Maple

September 13 2019

15 18

I’m very pleased to announce that we have just released the Maple Companion mobile app for iOS and Android phones. As its name implies, this free app is a complement to Maple. You can use it to take pictures of math you find out in the wild (e.g. in your handwritten notes, on a blackboard, in a textbook), and bring that math into Maple so you can get to work.

The Maple Companion lets you:

Avoid the mistakes that can occur when transcribing mathematical expressions into Maple manually
Save time when entering multiple equations into Maple, such as when you are checking your homework or pulling information from a reference book
Push math you’ll need later into Maple now, even if you don’t have your computer handy

The Maple Companion is an idea we started playing with recently. We believe it has interesting potential as a tool to help students learn math, and we’d really like your feedback to help shape its future direction. This first release is a step towards that goal, so you can try it out and start thinking about what else you would like to see from an app like this. Should it bring in entire documents? Integrate with tutors and Math Apps? Help students figure out where they went wrong when solving a problem? Let us know what you think!

Visit Maple Companion to learn more, link to the app stores so you can download the app, and access the feedback form. And of course, you are also welcome to give us your ideas in the comment section of this post.

Remove Sections in Maple Document

by: Neil K 40 Product: Maple

August 22 2019

8 3

I faced the issue of having to remove sections from a maple document in order to export to a pdf without the indentation and lines that come when you export documents with sections. Here is small tool I wrote that removes all sections in a maple document. It takes a target file as the first argument and writes that file without sections to the destination file specified as the second argument.

>

RemoveSection := module()
    local ModuleApply := proc(target, destination)
        XMLTools:-WriteFile(destination, subsindets(XMLTools:-ReadFile(target), ':-specfunc'(_XML_Section), section_handler));
    end proc;

    local section_handler := proc(s)
        local partresult := remove(type,[op(s)],`=`);
        return op(subsindets(partresult, ':-specfunc'(_XML_Title), f -> `_XML_Presentation-Block`("",_XML_Group("view"="presentation","inline-output"="false",
            "applyint"="true","applyrational"="true","applyexponent"="false","",_XML_Input(op(f))))));
    end proc;
end module:

>	#RemoveSection takes two arguments the first is the target file and the second is the destination file where the target file will be written without sections

>

Download RemoveSection.mw

An improved approximation of the Inverse Normal...

by: mmcdara 7891 Product: Maple

August 15 2019

3 20

Seeking for fast approximate formulas to compute (a huge number of) quantiles of a Gaussian random variable (here the standard one, but its extension to any Gaussian RV is straightforward), I found a few of them in the Abramowitz and Stegun book, page 933, relations 26.2.22 and 26.2.23.
Each approximation model is expressed as a rational fraction, the second one being the more accurate.
Each model depends on (respectively 4 and 6) parameters that are estimated (I guess it was done this way) through a least-square-like method.

See here for an online access http://people.math.sfu.ca/~cbm/aands/page_933.htm.

These approximation, and specially the most accurate one (formula 26.2.23) seem to be still widely used today⁽¹⁾ (see for instance https://www.johndcook.com/blog/normal_cdf_inverse/ ).

As an amusement I decided to compute the best fit by using the Statistics:-NonLinearFit procedure and a sample of (probability, quantile) points where probability ranges in [0.5, 1-1/1000] (the range used in formulas 26.2.22 and 26.2.23 is (0, 0.5] but this is not a point).
Surprisingly Statistics:-NonLinearFit returned, for the two formulas, parameter estimations substantially different from the one given in the Abramowitz & Stegun's book. A reason could be that the points I used when I did the fits weren't the one they used (unfortunately they give no informations about this).

More interesting, whatever the formula I refitted, NonLinearFit produced an approximation whose the absolute error was smaller by about two orders of magnitude to the onesprovided by Abramowitz and Stegun.
For instance they wrote that the most accurate formula (26.2.23) had an absolute approximation error less than 4.5*10^-4as I obtained a value around 10^-6!

(1) To get an idea of the persistence of the use of the formula 26.2.23, just type the value 2.515517 of its parameter c[0] in any search engine.

In the plots below the gray rectangle refers to the region where the approximate ICDF is used for extrapolation.

>

restart:

>	with(Statistics):

>	cdf := unapply(evalf(CDF(Normal(0, 1), x)), x): X := [seq(0..5, 0.1)]: A := cdf~(X): T := alpha -> sqrt(-2*log(1-alpha)): q := Quantile~(Normal(0, 1), A): Aq := convert([A,q], Matrix)^+:

>

r := 1:

J := z -> z - add(a__||k*z^k, k=0..r)/(1+add(b__||k*z^k, k=1..r+1)):

model := J(T(alpha)):
NL_fit := unapply(NonlinearFit(model, Aq, alpha), alpha);

# these lines are for estimating the performances
B := Sample(Uniform(0.5, 1), 10^4):
CodeTools:-Usage(Quantile~(Normal(0, 1), B)):
CodeTools:-Usage(Quantile~(Normal(0, 1), B, numeric)):
CodeTools:-Usage(NL_fit~(B)):
#-----------------------------------------------------
Y := [seq(0..6, 0.01)]:
B := cdf~(Y):
R1 := Quantile~(Normal(0, 1), B, numeric):
R2 := NL_fit~(B):

plots:-display(
ScatterPlot(R1, log[10]~(abs~(R2-~R1)), legend=mmcdara, color=red, gridlines=true, size=[700, 400]),
plottools:-rectangle([max(X), log[10]~(min(abs~(R2-~R1)))], [max(Y), log[10]~(max(abs~(R2-~R1)))], color=gray, transparency=0.6)
);

proc (alpha) options operator, arrow; (-2*ln(1-alpha))^(1/2)-(HFloat(2.5454311687345044)+HFloat(0.8058592540791468)*(-2*ln(1-alpha))^(1/2))/(1+HFloat(1.4689746699940707)*(-2*ln(1-alpha))^(1/2)-HFloat(0.34455942407858625)*ln(1-alpha)) end proc

memory used=170.31MiB, alloc change=76.01MiB, cpu time=3.06s, real time=3.05s, gc time=54.87ms

memory used=171.59MiB, alloc change=256.00MiB, cpu time=3.12s, real time=3.03s, gc time=154.77ms

memory used=8.24MiB, alloc change=0 bytes, cpu time=95.00ms, real time=95.00ms, gc time=0ns

>

r := 2:

J := z -> z - add(a__||k*z^k, k=0..r)/(1+add(b__||k*z^k, k=1..r+1)):

model := J(T(alpha)):
NL_fit := unapply(NonlinearFit(model, Aq, alpha), alpha);

# these lines are for estimating the performances
B := Sample(Uniform(0.5, 1), 10^4):
CodeTools:-Usage(Quantile~(Normal(0, 1), B)):
CodeTools:-Usage(Quantile~(Normal(0, 1), B, numeric)):
CodeTools:-Usage(NL_fit~(B)):
#-----------------------------------------------------

Y := [seq(0..6, 0.01)]:
B := cdf~(Y):
R1 := Quantile~(Normal(0, 1), B, numeric):
R2 := NL_fit~(B):

plots:-display(
ScatterPlot(R1, log[10]~(abs~(R2-~R1)), legend=mmcdara, color=red, gridlines=true, size=[700, 400]),
plottools:-rectangle([max(X), log[10]~(min(abs~(R2-~R1)))], [max(Y), log[10]~(max(abs~(R2-~R1)))], color=gray, transparency=0.6)
);

proc (alpha) options operator, arrow; (-2*ln(1-alpha))^(1/2)-(HFloat(2.9637294443959394)+HFloat(4.527738737327481)*(-2*ln(1-alpha))^(1/2)-HFloat(0.9571637188191973)*ln(1-alpha))/(1+HFloat(3.472400103322335)*(-2*ln(1-alpha))^(1/2)-HFloat(3.426536241250657)*ln(1-alpha)+HFloat(0.08875278252087411)*(-2*ln(1-alpha))^(3/2)) end proc

memory used=170.09MiB, alloc change=32.00MiB, cpu time=3.29s, real time=3.11s, gc time=268.60ms

memory used=170.85MiB, alloc change=0 bytes, cpu time=3.23s, real time=3.10s, gc time=201.52ms
memory used=10.76MiB, alloc change=0 bytes, cpu time=127.00ms, real time=127.00ms, gc time=0ns

>

# Optimized "r=2" computation

z_fit := simplify(subs(alpha=-exp(-(1/2)*z^2)+1, NL_fit(alpha))) assuming z > 0:
z_fit := unapply(convert~(%, horner), z);

p := proc(alpha)
  local z:
  z := sqrt(-2*log(1-alpha)):
  z_fit(z):
end proc:

R3 := CodeTools:-Usage(p~(B)):

plots:-display(
  ScatterPlot(R1, log[10]~(abs~(R2-~R1)), legend=mmcdara, color=red, gridlines=true, size=[700, 400]),
  plottools:-rectangle([max(X), log[10]~(min(abs~(R2-~R1)))], [max(Y), log[10]~(max(abs~(R2-~R1)))], color=gray, transparency=0.6)
);

memory used=1.67MiB, alloc change=0 bytes, cpu time=14.00ms, real time=15.00ms, gc time=0ns

AS stands for Abramowith & Stegun

>

J_AS := unapply(normal(eval(J(t), [a__0=2.515517, a__1=0.802853, a__2=0.010328, b__1=1.432788, b__2=0.189269, b__3=0.001308])), t):
J_AS(t);

# for comparison:

print():
z_fit := simplify(subs(alpha=-exp(-(1/2)*z^2)+1, NL_fit(alpha))) assuming z > 0:
map(sort, %, z);

plot([z_fit(z), J_AS(z)], z=0.5..1, color=[blue, red], legend=[mmcdara, Abramowitz_Stegun], gridlines=true);

print():
R2_AS := CodeTools:-Usage(J_AS~(T~(B))):
print():

plots:-display(
  ScatterPlot(R1, log[10]~(abs~(R2_AS-~R1)), legend=Abramowitz_Stegun, gridlines=true, size=[700, 400]),
  ScatterPlot(R1, log[10]~(abs~(R2-~R1)), legend=mmcdara, color=red),
  plottools:-rectangle([max(X), log[10]~(min(abs~(R2-~R1)))], [max(Y), log[10]~(max(abs~(R2-~R1)))], color=gray, transparency=0.6)
);

(0.1308000000e-2*t^4+.1892690000*t^3+1.422460000*t^2+.1971470000*t-2.515517000)/(0.1308000000e-2*t^3+.1892690000*t^2+1.432788000*t+1.)

(0.8875278252e-1*z^4+1.713268121*z^3+2.993818244*z^2-3.527738737*z-2.963729444)/(0.8875278252e-1*z^3+1.713268121*z^2+3.472400103*z+1.)

memory used=2.92MiB, alloc change=0 bytes, cpu time=25.00ms, real time=25.00ms, gc time=0ns

>

Download InverseNormalCDF.mw

Preparing a parametrized surface for 3D printing...

by: Rouben Rostamian 8506 Product: Maple

August 15 2019

8 4

Although the graph of a parametrized surface can be viewed and manipulated on the computer screen as a surface in 3D, it is not quite suitable for printing on a 3D printer since such a surface has zero thickness, and thus it does not correspond to physical object.

To produce a 3D printout of a surface, it needs to be endowed with some "thickness". To do that, we move every point from the surface in the direction of that point's nomral vector by the amount ±T/2, where T is the desired thickness. The locus of the points thus obtained forms a thin shell of thickness T around the original surface, thus making it into a proper solid. The result then may be saved into a file in the STL format and be sent to a 3D printner for reproduction.

The worksheet attached to this post provides a facility for translating a parametrized surface into an STL file. It also provides a command for viewing the thickened object on the screen. The details are documented within that worksheet.

Here are a few samples. Each sample is shown twice—one as it appears within Maple, and another as viewed by loading the STL file into MeshLab which is a free mesh viewing/manipulation software.

Here is the worksheet that produced these: thicken.mw

Structural Dynamics with Maplesim

by: Lenin Araujo Castillo 1790 Product: MapleSim 2018

August 10 2019

4 0

Analysis in Dynamics of Structures with Maplesim for Engineering
Here is the power of Maplesim in modeling and simulation. With Maplesim you can model structures at rest and dynamics. Considering real patterns of our world for better optimization.Project developed for students of Civil Engineering, Architecture, Mechatronics and all those professional careers related to structures.

CIMAC_UNALM_2019.pdf

Lenin Araujo Castillo

Ambassador of Maple

Serpentine paths in matrices and generation of...

by: Kitonum 21435 Product: Maple 2018

August 07 2019

3 0

This post is closely related to the previous one https://www.mapleprimes.com/posts/210930-Numbrix-Puzzle-By-The-Branch-And-Bound-Method which presents the procedure NumbrixPuzzle that allows you to effectively solve these puzzles (the text of this procedure is also available in the worksheet below).
This post is about generating these puzzles. To do this, we need the procedure SerpentinePaths (see below) , which allows us to generate a large number of serpentine paths in a matrix of a specified size, starting with a specified matrix element. Note that for a square matrix of the order n , the number of such paths starting from [1,1] - position is the sequence https://oeis.org/search?q=1%2C2%2C8%2C52%2C824&language=english&go=Search .

The required parameter of SerpentinePaths procedure is the list S , which defines the dimensions of the matrix. The optional parameter is the list P - this is the position of the number 1 (by default P=[1,1] ).
As an example below, we generate 20 puzzles of size 6 by 6. In exactly the same way, we can generate the desired number of puzzles for matrices of other sizes.

>

restart;

>

SerpentinePaths:=proc(S::list, P::list:=[1,1])
local OneStep, A, m, F, B, T, a;

OneStep:=proc(A::listlist)
local s, L, B, T, k, l;

s:=max[index](A);
L:=[[s[1]-1,s[2]],[s[1]+1,s[2]],[s[1],s[2]-1],[s[1],s[2]+1]];
T:=table(); k:=0;
for l in L do
if l[1]>=1 and l[1]<=S[1] and l[2]>=1 and l[2]<=S[2] and A[op(l)]=0 then k:=k+1; B:=subsop(l=a+1,A);
T[k]:=B fi;
od;
convert(T, list);
end proc;
A:=convert(Matrix(S[], {(P[])=1}), listlist);
m:=S[1]*S[2]-1;
B:=[$ 1..m];
F:=LM->ListTools:-FlattenOnce(map(OneStep, `if`(nops(LM)<=30000,LM,LM[-30000..-1])));
T:=[A];
for a in B do
T:=F(T);
od;
map(convert, T, Matrix);

end proc:

>

NumbrixPuzzle:=proc(A::Matrix)
local A1, L, N, S, MS, OneStepLeft, OneStepRight, F1, F2, m, L1, p, q, a, b, T, k, s1, s, H, n, L2, i, j, i1, j1, R;
uses ListTools;
S:=upperbound(A); N:=nops(op(A)[3]); MS:=`*`(S);
A1:=convert(A, listlist);
for i from 1 to S[1] do
for j from 1 to S[2] do
for i1 from i to S[1] do
for j1 from 1 to S[2] do
if A1[i,j]<>0 and A1[i1,j1]<>0 and abs(A1[i,j]-A1[i1,j1])<abs(i-i1)+abs(j-j1) then return `no solutions` fi;
od; od; od; od;
L:=sort(select(e->e<>0, Flatten(A1)));
L1:=[`if`(L[1]>1,seq(L[1]-k, k=0..L[1]-2),NULL)];
L2:=[seq(seq(`if`(L[i+1]-L[i]>1,L[i]+k,NULL),k=0..L[i+1]-L[i]-2), i=1..nops(L)-1), `if`(L[-1]<MS,seq(L[-1]+k,k=0..MS-L[-1]-1),NULL)];
OneStepLeft:=proc(A1::listlist)
local s, M, m, k, T;
uses ListTools;
s:=Search(a, Matrix(A1));
M:=[[s[1]-1,s[2]],[s[1]+1,s[2]],[s[1],s[2]-1],[s[1],s[2]+1]];
T:=table(); k:=0;
for m in M do
if m[1]>=1 and m[1]<=S[1] and m[2]>=1 and m[2]<=S[2] and A1[op(m)]=0 then k:=k+1; T[k]:=subsop(m=a-1,A1);
fi;
od;
convert(T, list);
end proc;
OneStepRight:=proc(A1::listlist)
local s, M, m, k, T, s1;
uses ListTools;
s:=Search(a, Matrix(A1)); s1:=Search(a+2, Matrix(A1));
M:=[[s[1]-1,s[2]],[s[1]+1,s[2]],[s[1],s[2]-1],[s[1],s[2]+1]];
T:=table(); k:=0;
for m in M do
if m[1]>=1 and m[1]<=S[1] and m[2]>=1 and m[2]<=S[2] and A1[op(m)]=0 and `if`(a+2 in L, `if`(is(abs(s1[1]-m[1])+abs(s1[2]-m[2])>1),false,true),true) then k:=k+1; T[k]:=subsop(m=a+1,A1);
fi;
od;
convert(T, list);
end proc;
F1:=LM->ListTools:-FlattenOnce(map(OneStepLeft, LM));
F2:=LM->ListTools:-FlattenOnce(map(OneStepRight, LM));
T:=[A1];
for a in L1 do
T:=F1(T);
od;
for a in L2 do
T:=F2(T);
od;
R:=map(t->convert(t,Matrix), T);
if nops(R)=0 then return `no solutions` else R fi;
end proc:

Simple examples

>

SerpentinePaths([3,3]); # All the serpentine paths for the matrix 3x3, starting with [1,1]-position
SerpentinePaths([3,3],[1,2]); # No solutions if the start with [1,2]-position
SerpentinePaths([4,4]): # All the serpentine paths for the matrix 4x4, starting with [1,1]-position
nops(%);
nops(SerpentinePaths([4,4],[1,2])); # The number of all the serpentine paths for the matrix 4x4, starting with [1,2]-position
nops(SerpentinePaths([4,4],[2,2])); # The number of all the serpentine paths for the matrix 4x4, starting with [2,2]-position

(1)

Below we find 12,440 serpentine paths in the matrix 6x6 starting from various positions (the set L )

>	k:=0: n:=6: for i from 1 to n do for j from i to n do k:=k+1; S[k]:=SerpentinePaths([n,n],[i,j])[]; od: od: L1:={seq(S[i][], i=1..k)}: L2:=map(A->A^%T, L1): L:=L1 union L2: nops(L);

(2)

Further, using the list L, we generate 20 examples of Numbrix puzzles with the unique solutions

>

T:='T':
N:=20:
M:=[seq(L[i], i=combinat:-randcomb(nops(L),N))]:
for i from 1 to N do
for k from floor(n^2/4) do
T[i]:=Matrix(n,{seq(op(M[i])[3][j], j=combinat:-randcomb(n^2,k))});
if nops(NumbrixPuzzle(T[i]))=1 then break; fi;
od: od:
T:=convert(T,list):
T1:=[seq([seq(T[i+j],i=1..5)],j=[0,5,10,15])]:
DocumentTools:-Tabulate(Matrix(4,5, (i,j)->T1[i,j]), fillcolor = "LightYellow", width=95):

The solutions of these puzzles

>	DocumentTools:-Tabulate(Matrix(4,5, (i,j)->NumbrixPuzzle(T1[i,j])[]), fillcolor = "LightYellow", width=95):

>

For some reason, these 20 examples and their solutions did not load here.

Edit. I separately inserted these generated 20 puzzles as a picture:

Download SerpPathsinMatrix.mw

Membership Anniversary For Ten Years (10)

by: Fereydoon_Shekofte 223 Product: MaplePrimes

August 07 2019

1 0

One decade on MaplePrimes

Always I Wish A Wonderful Year In Mathematics As Annus mirabilis

Fereydoon Shekofte

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

MaplePrimes Posts

MaplePrimes Posts are for sharing your experiences, techniques and opinions about Maple, MapleSim and related products, as well as general interests in math and computing.

Latest Post
Latest Posts Feed

Stiffness matrix implementation of a finite element...

Free vibration of a cracked composite beam

Splitting PDE parameterized symmetries and parameter...

Pole simulation with tensions in MapleSim

Maple May benefit from starting a Stack Exchange...

Integral Transforms (revamped) and PDE

Graph Theory and Pokémon Application

Arabic Cipher

Introducing the Maple Companion app!

Remove Sections in Maple Document

An improved approximation of the Inverse Normal...

Preparing a parametrized surface for 3D printing...

Structural Dynamics with Maplesim

Serpentine paths in matrices and generation of...

Membership Anniversary For Ten Years (10)

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again

Ask a Question

Create a Post

MaplePrimes Posts

MaplePrimes Posts are for sharing your experiences, techniques and opinions about Maple, MapleSim and related products, as well as general interests in math and computing.

Latest Post Latest Posts Feed

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Latest Post
Latest Posts Feed