henrylyl

20 Reputation

2 Badges

11 years, 36 days

MaplePrimes Activity


These are questions asked by henrylyl

Does anyone has any maple worksheet that generate surface using the PDE method described in this article?  I am trying to learn this method but I am not familiar with the mathematics to do it although the paper gives some description of it.  I hope someone can demonstrate the procedure in Maple.  Thanks

http://www.researchgate.net/publication/259095177_Automatic_shape_optimisation_of_pharmaceutical_tablets_using_Partial_Differential_Equations/file/72e7e52a87ed129d4a.pdf

Please see the attached worksheet for the error at the end, any idea what happens?

Besides, anyone know how to make plot array with each element has a plot and scatterplot in it?  I made my plot by making table first and then type the cmd individually.

Raw Data

restart

avicel102_aero200 := Matrix(12, 3, {(1, 1) = 50, (1, 2) = 1.6621, (1, 3) = .22, (2, 1) = 100, (2, 2) = 3.33221, (2, 3) = .103, (3, 1) = 150, (3, 2) = 4.95533, (3, 3) = 0.43e-1, (4, 1) = 200, (4, 2) = 5.62147, (4, 3) = 0.37e-1, (5, 1) = 50, (5, 2) = 1.88627, (5, 3) = .207, (6, 1) = 50, (6, 2) = 1.90375, (6, 3) = .199, (7, 1) = 100, (7, 2) = 4.04708, (7, 3) = 0.89e-1, (8, 1) = 100, (8, 2) = 4.20413, (8, 3) = 0.87e-1, (9, 1) = 150, (9, 2) = 5.50509, (9, 3) = 0.41e-1, (10, 1) = 150, (10, 2) = 4.90453, (10, 3) = 0.49e-1, (11, 1) = 200, (11, 2) = 6.22001, (11, 3) = 0.16e-1, (12, 1) = 200, (12, 2) = 6.24435, (12, 3) = 0.23e-1}); avicel102_syloid := Matrix(13, 3, {(1, 1) = 50, (1, 2) = 1.76578, (1, 3) = .235, (2, 1) = 100, (2, 2) = 4.08314, (2, 3) = .128, (3, 1) = 150, (3, 2) = 4.85921, (3, 3) = 0.73e-1, (4, 1) = 200, (4, 2) = 5.66786, (4, 3) = 0.51e-1, (5, 1) = 50, (5, 2) = 2.06221, (5, 3) = .224, (6, 1) = 50, (6, 2) = 2.17294, (6, 3) = .216, (7, 1) = 100, (7, 2) = 4.03728, (7, 3) = .109, (8, 1) = 100, (8, 2) = 3.67978, (8, 3) = .116, (9, 1) = 150, (9, 2) = 5.78427, (9, 3) = 0.6e-1, (10, 1) = 150, (10, 2) = 5.43599, (10, 3) = 0.59e-1, (11, 1) = 200, (11, 2) = 6.26332, (11, 3) = 0.38e-1, (12, 1) = 200, (12, 2) = 6.10814, (12, 3) = 0.39e-1, (13, 1) = 200, (13, 2) = 6.83857, (13, 3) = 0.25e-1}); avicel301_101_aero200 := Matrix(12, 3, {(1, 1) = 50, (1, 2) = 1.106, (1, 3) = .2597543205, (2, 1) = 50, (2, 2) = 1.06, (2, 3) = .2564849561, (3, 1) = 50, (3, 2) = 1.129, (3, 3) = .2552025392, (4, 1) = 100, (4, 2) = 2.865, (4, 3) = .1280363098, (5, 1) = 100, (5, 2) = 2.793, (5, 3) = .127659333, (6, 1) = 100, (6, 2) = 3.088, (6, 3) = .1347147061, (7, 1) = 150, (7, 2) = 4.216, (7, 3) = 0.802489497e-1, (8, 1) = 150, (8, 2) = 4.044, (8, 3) = 0.838493001e-1, (9, 1) = 150, (9, 2) = 4.071, (9, 3) = 0.855901289e-1, (10, 1) = 200, (10, 2) = 4.826, (10, 3) = 0.627298821e-1, (11, 1) = 200, (11, 2) = 4.924, (11, 3) = 0.614759445e-1, (12, 1) = 200, (12, 2) = 4.76, (12, 3) = 0.61001365e-1})

NULL

 

with(LinearAlgebra):

X[1] := Column(avicel102_aero200, 3):

X[2] := Column(avicel102_syloid, 3):

X[3] := Column(avicel301_101_aero200, 3):

S := Array(1 .. 3):

S[1] := ScatterPlot(X[1], Y[1], labels = [Porocity, Tensile*Strength], title = Compactability*of*Avicel102*with*Aerosil200):

with(plots):NULL

 

vSST := Array(1 .. 3):

tpfit := proc (x, y) options operator, arrow; Fit(a*exp(-k*p)+c, x, y, p, output = [leastsquaresfunction, residualsumofsquares, degreesoffreedom]) end proc:SST := proc (y) options operator, arrow; Variance(y)*(NumElems(y)-1) end proc:NULL

for i to 3 do fnTP[i] := tpfit(X[i], Y[i])[1] end do:

for i to 3 do fnPC[i] := Fit(c*p^2+b*p+a, P[i], X[i], p, output = [leastsquaresfunction, residualsumofsquares, degreesoffreedom])[1] end do:

NULL

for i to 3 do vSST[i] := SST(Y[i]) end do:

for i to 3 do pSST[i] := SST(X[i]) end do:

NULL

for i to 3 do vSSR[i] := tpfit(X[i], Y[i])[2] end do:

for i to 3 do Rsquare[i] := 1-vSSR[i]/vSST[i] end do:

NULL

``

s := Array(1 .. 3):

for i to 3 do s[i] := plot(fnTP[i], p = 0 .. .26, legend = typeset(fnTP[i], R^2, "= ", Rsquare[i])) end do

for i to 3 do ps[i] := plot(fnPC[i], p = 0 .. 200, legend = typeset(fnPC[i], R^2, "= ", pRsquare[i])) end do

G := Array(1 .. 3):

G[1] := [S[1], s[1]]:

pG := [[PS[1], ps[1]], [PS[2], ps[2]], [PS[3], ps[3]]]:NULL

NULL

Table 1: Compactability Curve

display(G[1])

 

NULL

display(G[2])

 

NULL

display(G[3])

 

NULL

NULL

NULL

Table 2: Compressibility Curve

display(pG[1])

 

NULL

display(pG[2])

 

NULL

display(pG[3])

 

NULL

``

ff := proc (x) options operator, arrow; unapply(x, p) end proc:

gtp := map(ff, fnTP):

gpc := map(ff, fnPC):

TS1 := proc (pressure) options operator, arrow; (`@`(gtp[1], gpc[1]))(pressure/Unit('MPa'))*Unit('MPa') end proc:

TS2 := proc (pressure) options operator, arrow; (`@`(gtp[2], gpc[2]))(pressure/Unit('MPa'))*Unit('MPa') end proc:

TS3 := proc (pressure) options operator, arrow; (`@`(gtp[3], gpc[2]))(pressure/Unit('MPa'))*Unit('MPa') end proc:

por1 := proc (p) options operator, arrow; gpc[1](p/Unit('MPa')) end proc:

(1)

por1(x*Unit('MPa'))

HFloat(0.36333333333333295)-HFloat(0.0036139999999999957)*x+HFloat(9.666666666666651e-6)*x^2

(2)

por1(x);

HFloat(0.36333333333333295)-HFloat(0.0036139999999999957)*x/Units:-Unit('MPa')+HFloat(9.666666666666651e-6)*x^2/Units:-Unit('MPa')^2

(3)

``

NULL

 

(4)

ts := 2/3*(10*P/(Pi*D^2*(2.84*t/D-.126*t/W+3.15*W/D+0.1e-1))):

Target Weight of Tablet, m__t:

m__t := 700*Unit('mg')

700*Units:-Unit('mg')

(5)

Tablet True Density:NULL

`ρ__true1` := 1.8*Unit('g'/'cm'^3)

1800.0*Units:-Unit(('kg')/('m')^3)

(6)

Width of tablet w;

d := 7.747*Unit('mm')

7.747*Units:-Unit('mm')

(7)

Length of tablet, l:

l := 17.018*Unit('mm')

17.018*Units:-Unit('mm')

(8)

Cup Depth, h:

hh := 1.6764*Unit('mm')

1.6764*Units:-Unit('mm')

(9)

Cup Volume", `v__cup`"

v__cup := 119.789438*Unit('mm'^3)

119.789438*Units:-Unit(('mm')^3)

(10)

Die cross sectional area, v__die:

a__die := 118.9610524*Unit('mm'^2)

118.9610524*Units:-Unit(('mm')^2)

(11)

Equation to solve tablet thickness using tablet dimension and porocity:

Tablet bulk density:

`ρ__bulk` := `ρ__true`*sf

`ρ__true`*sf

(12)

Tablet volume:``

V__tablet := 2*V__cup+A__die*(t-2*H)

2*V__cup+A__die*(t-2*H)

(13)

NULL

Tablet bulk Density:``

`ρ__bulk` := m__tablet/V__tablet

m__tablet/(2*V__cup+A__die*(t-2*H))

(14)

NULL

eq1 := m__tablet/(2*V__cup+A__die*(t-2*H)) = `ρ__true`*sf

m__tablet/(2*V__cup+A__die*(t-2*H)) = `ρ__true`*sf

(15)

Function for thickness:

thickness := unapply(simplify(solve(eq1, t), size), sf, A__die, H, V__cup, `ρ__true`, m__tablet)

proc (sf, A__die, H, V__cup, `ρ__true`, m__tablet) options operator, arrow; (2*sf*(A__die*H-V__cup)*`ρ__true`+m__tablet)/(A__die*sf*`ρ__true`) end proc

(16)

simplify(thickness(1-por[1](100*Unit(MPa)), a__die, hh, v__cup, `ρ__true1`, m__t))

HFloat(0.005719286791216391)*Units:-Unit('m')

(17)

t2 := proc (cc, A__die, H, V__cup, `ρ__true`, m__tablet) options operator, arrow; simplify(thickness(1-cc, A__die, H, V__cup, `ρ__true`, m__tablet)) end proc:

``

``

``

Tablet tensile strength:NULL

eq2 := `σ__t` = ts

`σ__t` = (20/3)*P/(Pi*D^2*(2.84*t/D-.126*t/W+3.15*W/D+0.1e-1))

(18)

Break force of tablet (N)

breakforce := unapply(subs(W = t-2*h, solve(eq2, P)), `σ__t`, D, h, t)

proc (`σ__t`, D, h, t) options operator, arrow; 0.9424777961e-3*`σ__t`*D*(5.*D*(t-2*h)-63.*t*D+1575.*(t-2*h)^2+1420.*t*(t-2*h))/(t-2*h) end proc

(19)

 

kp := simplify(breakforce(TS[1](100*Unit('MPa')), d, hh, 0.5719e-2*Unit('m')))

HFloat(286.99924404612904)*Units:-Unit('N')

(20)

convert(kp, units, kilopond)

HFloat(29.265778226624693)*Units:-Unit('kilopond')

(21)

t2(por[1](100*Unit('MPa')), a__die, hh, v__cup, `ρ__true1`, m__t)

HFloat(0.005719286791216391)*Units:-Unit('m')

(22)

``

hh

1.6764*Units:-Unit('mm')

(23)

simplify(breakforce(TS[1](100*Unit('MPa')), d, hh, t2(por[1](100*Unit('MPa')), a__die, hh, v__cup, `ρ__true1`, m__t)))

HFloat(287.02451456724083)*Units:-Unit('N')

(24)

Encounter problem while chaning variable

 BF:= 
proc (`σ__t`, D, p, A__die, hhh, V__cup, `ρ__true`, m__tablet) options operator, arrow; breakforce(`σ__f`, D, hhh, t2(p, A__die, hhh, V__cup, `ρ__true`, m__tablet)) end proc:

BF(TS1(100*Unit('MPa')), d, por1(100*Unit('MPa')), a__die, hh, v__cup, `ρ__true1`, m__t)

BF(HFloat(3.6712504621161846)*Units:-Unit('MPa'), 7.747*Units:-Unit('mm'), HFloat(0.09859999999999991), 118.9610524*Units:-Unit(('mm')^2), 1.6764*Units:-Unit('mm'), 119.789438*Units:-Unit(('mm')^3), 1420.00*Units:-Unit(('kg')/('m')^3), 667*Units:-Unit('mg'))

(25)

evalf(TS[1](100*Unit(MPa)))

HFloat(3.6712504621161846)*Units:-Unit('MPa')

(26)

plot(breakforce(TS1(x), d, hh, t2(por1(x), a__die, hh, v__cup, `ρ__true1`, m__t)), x = Unit('MPa') .. 200*Unit('MPa'), useunits = true)

Error, (in Units:-Standard:-+) the units `1` and `1/MPa` have incompatible dimensions

 

plot(t2(por1(x), a__die, hh, v__cup, `ρ__true1`, m__t), x = Unit('MPa') .. 200*Unit('MPa'), useunits = true)

Error, (in Units:-Standard:-+) the units `1` and `1/MPa` have incompatible dimensions

 

plot(evalf(por1(xx)), xx = 10 .. 200)

 

but por1 works fine on its own with units

por1(Units:-Standard:-`*`(200, Unit('MPa')))

HFloat(0.02719999999999989)

(27)

plot(evalf(por1(xx)), xx = 10 .. 200);

 

Unloading Units:-Standard

``

``

``

``

``

``

``

``

``

``

``

``

``

``

 

Fitting_worksheet.mw

I also used a lot for for loop in this worksheet, can anyone suggest another way of doing it? And I think I initiated many Array(1..3) that is of no use.  But if I do not initiate and do the following

A[1]:=a : 

A[2]:=b:

A[3]:=c:

A;

output: A

Only if I initiate A:=Array(1..3) and then define the element I can see the complete A array at the end.

 

I am trying to write a procedure that have the following input:  expr, xdata,ydata. [dependent variables list]

I want a procedure to calculate Rsquare following the method described in 

http://www.mapleprimes.com/questions/120014-How-To-Calculate-Coefficient-Of-Determination

 

But I want the proc to take any abitary function and be able to identify the  dependent variable...how can I approach it? what function can I use to match [dependent variable list] in expr??

 

Thanks

 

 

I have a question regarding maple reader.

 

It seems like in maple reader, one cannot interact with the output of the Explore function.... which mean I am forced to make my own component if I want to creat application for my co-worker to use..this is some what annoying.  Is there another way around?

 

please open the attached file in maple player and you will see the output created by Explore gives error message.  mapleplayer_test.mw

 

Download mapleplayer_test.mw

.  

Hi Mapleprimes community,

 

I am a maple novice and was trying to write a procedure but cannot tell what my syntax error is.  I would greatly appreciate if anyone can point out what I did wrong or any improvement in the code that I can make.  Thanks

SLRrepeatedsample_pr.mw

 

Henry

Page 1 of 1