mmcdara

If I'm not mistaken...

Answered: mmcdara 7891

May 28 2023

0 0

Based upon Descarte's tule of signs, and assuming gamma > 0, I would say there is only one positive root.
If gamma < 0 I can' be that conclusive.
Only_1_positive_root.mw

A solution...

Answered: mmcdara 7891

May 28 2023

1 3

Here is a solution (an example about how to place the legend within the plot and to build a composite legend made of line style and symbol ifs proposed in the file):
A_solution.mw

Download A_solution.mw

You can also refer to the answer I gave to your earlier question question https://www.mapleprimes.com/questions/236446-How-To-To-Draw-The-Horizontal-And-Veritcal

Here are two ways:...

Answered: mmcdara 7891

May 27 2023

1 6

The first option is immediate, the second one places t "legend-ike box" within the plot window itself.

>

restart:

>	N := (-18sqrt(2)(M-4/3)(x^2+2)^2arctan((1/2)xsqrt(2))-9Pi(M-4/3)(x^2+2)^2sqrt(2)+(-36M+48)x^3+(-120M+96)x)/(4*(x^2+2)^2)

(1/4)*(-18*2^(1/2)*(M-4/3)*(x^2+2)^2*arctan((1/2)*x*2^(1/2))-9*Pi*(M-4/3)*(x^2+2)^2*2^(1/2)+(-36*M+48)*x^3+(-120*M+96)*x)/(x^2+2)^2

(1)

>

p := plot(
  [subs(M = 1.3015, N)]
  , x = -5 .. 5
  , labels = ["x", "N"]
  , axes = boxed
  , labeldirections = ["horizontal", "horizontal"]
  , colour = red
  , linestyle = solid
  , legend = ["N= 1, M= 1.3015"]
  , legendstyle = [location = bottom, font = ["TIMES", "italic", 15]]
):

ymin, ymax := (min, max)(plottools:-getdata(p)[3][.., 2]):
eps   := 0.1:
Vline := plot([[0, ymin-eps], [0, ymax+eps]], color=black):
Hline := plot([[-5, 0], [5, 0]], color=black):

plots:-display(p, Vline, Hline);

>

p := plot(
  [subs(M = 1.3015, N)]
  , x = -5 .. 5
  , labels = ["x", "N"]
  , axes = boxed
  , labeldirections = ["horizontal", "horizontal"]
  , colour = red
  , linestyle = solid
  , gridlines = false
):

ymin, ymax := (min, max)(plottools:-getdata(p)[3][.., 2]):
eps   := 0.1:
Vline := plot([[0, ymin-eps], [0, ymax+eps]], color=black):
Hline := plot([[-5, 0], [5, 0]], color=black):

Left := -4.8:
Right := -0.5:

Legend := plots:-display(
  plot(ymin, x=Left..Left+1, color=red)
  , plots:-textplot([Left+1, ymin, " N= 1, M= 1.3015"], align=right)
  , plottools:-rectangle([-4.9, ymin-eps], [-0.2, ymin+eps], color=white)
  #, plot([[Left-0.1, ymin-eps], [Right, ymin-eps], [Right, ymin+eps], [Left-0.1, ymin+eps], [Left-0.1, ymin-eps]], color=black)
):

plots:-display(p, Vline, Hline, Legend);

>

Download TwoOptions.mw

Waiting for concurrent proposals......

Answered: mmcdara 7891

May 27 2023

0 1

... here is a first answer.
The coding is simple and probably not that efficient vor (very) large matrices, Let's just say it's the start of the exchanges to come?

>

restart

>	with(LinearAlgebra):

>	M := RandomMatrix(6$2, generator=0..10)

(1)

>

# A rough code: it"s likely that adjustments must be done for large matrices
# in order to get an efficient code.

SmallestBelow := proc(M, withdiag)
  local coding, minval, here:
  coding := x -> piecewise(x = 0, 1, 0);
  map(x -> [lhs(x)], op(coding~(M))[3]);
  if withdiag then
    minval := min(seq(seq(M[i, j], j=1..i), i=1..numelems(M[1])));
    here   := select((x -> is(x[2] <= x[2])), map(x -> [lhs(x)], op(coding~(M))[3]));
  else
    minval := min(seq(seq(M[i, j], j=2..i), i=1..numelems(M[1])));
    here   := select((x -> is(x[2] < x[1])), map(x -> [lhs(x)], op(coding~(M))[3]));
  end if:
  minval, here
end proc:

# Strictly delow the diagonal

minval, here := SmallestBelow(M, false);

# Diagonal included

minval, here := SmallestBelow(M, true);

(2)

>

Download Min_value_and_locations.mw

Obvious solution:...

Answered: mmcdara 7891

May 26 2023

0 0

(x-4*y)^9+3*z reaches its maximum value when z=17 and when (x-4*y)^9 is maximum, thus when x=3 and y=0.
Maximum value = 3^9+3*17 = 19734

Didn't I already answered simi...

Answered: mmcdara 7891

May 25 2023

1 3

Didn't I already answered a similar question weeks ago ???

In case you didn't read my answer, which can be possible, here it is again.
In case you read it but didn't like it, I would have been happy to receive your comments instead of seeing it again under hidden form.
By the way, the code below works for any version of Maple.

>

restart;

>	#with(student):

>	U := a[0] + sum(-a[i]tanh(xi[n])^i, i = 1 .. 1) + sum(-b[i]tanh(xi[n])^(-i), i = 1 .. 1):

>	u(xi[n + 1]) := a[0] - a[1](tanh(xi[n]) + tanh(d))/(1 + tanh(xi[n])tanh(d)) - b[1](1 + tanh(xi[n])tanh(d))/(tanh(xi[n]) + tanh(d)):

>	u(xi[n - 1]) := a[0] - a[1](tanh(xi[n]) - tanh(d))/(1 - tanh(xi[n])tanh(d)) - b[1](1 - tanh(xi[n])tanh(d))/(tanh(xi[n]) - tanh(d)):

>	eq := c[1]^2diff(U, xi[n], xi[n])ln(1 + U/alpha) - beta(u(xi[n + 1]) - 2U - u(xi[n - 1])):

>	fin1 := simplify(numer(eq)):

>	fin := simplify(subs(tanh(xi[n]) = Psi, fin1)):

>

# As the degree(Psi*ln(Psi), Psi) returns FAIL, the idea is to replace ln(...) by
# something else, let's say H.
#
# Observe that EQ[1] is the identity 0=0

lnPsi := select(has, indets(fin, function), ln)[]:
fin_H := eval(fin, lnPsi=H):
deg   := degree(fin_H, Psi);
print():

for i from 0 to deg do
  collect(coeff(fin_H, Psi, i), H);
  EQ[i] := eval(%, H=lnPsi) = 0;
end do:

print~(EQ):

-2*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))*tanh(d)^2*b[1]*c[1]^2 = 0

(2*tanh(d)^4*b[1]*c[1]^2+2*tanh(d)^2*b[1]*c[1]^2+2*b[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))-2*tanh(d)^2*beta*b[1] = 0

(-2*tanh(d)^4*b[1]*c[1]^2+2*tanh(d)^2*a[1]*c[1]^2-2*tanh(d)^2*b[1]*c[1]^2-2*b[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))+2*tanh(d)^4*beta*b[1]-2*tanh(d)^2*beta*a[1]+2*beta*b[1] = 0

(-2*tanh(d)^4*a[1]*c[1]^2-2*tanh(d)^2*a[1]*c[1]^2+2*tanh(d)^2*b[1]*c[1]^2-2*a[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))+2*tanh(d)^4*beta*a[1]-2*tanh(d)^2*beta*b[1]+2*beta*a[1] = 0

(2*tanh(d)^4*a[1]*c[1]^2+2*tanh(d)^2*a[1]*c[1]^2+2*a[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))-2*tanh(d)^2*beta*a[1] = 0

-2*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))*tanh(d)^2*a[1]*c[1]^2 = 0

(1)

>

Download Degree_3_mmcdara.mw

@Kitonum's answer is likely the simp...

Answered: mmcdara 7891

May 23 2023

0 0

@Kitonum's answer is likely the simpler and the most efficient.

Nevertheless here is a slight variant:

for i from 3 to nops(n) do 
  if i::odd then
    f[i] := n[i-2]-3*n[i-1]+3*n[i]; 
  else
    f[i] := n[i-3]+4*n[i-2]-6*n[i-1]+n[i]; 
  end if:
end do:

@Rouben Rostamian Hi Rouben, Not...

Answered: mmcdara 7891

May 22 2023

2 2

@Rouben Rostamian

Hi Rouben,
Not in my intention to interfere in your exchange with the OP, but I suggest you to do something like this, whith your more recent version than my old 2015.2 (I was forced to rotate the plot to circumvent the lack of the orientation option in texplot)

>

restart;

>	f := theta -> 6000*sin(theta) - 3819.718635

(1)

>

pl1 := plots:-display(
    plot(f(t)+3819.718635, t=0..Pi, coords=polar, color=red, thickness=3)
  , seq(plot(r, t=0..Pi, coords=polar, color=gray, thickness=1), r in [seq](0..7000, 3819.718635/5))
  , plot(3819.718635, t=0..Pi, coords=polar, color=black, thickness=1)
  , seq(plot([[0, 0], [7000*cos(t*Pi), 7000*sin(t*Pi)]], color=gray), t in [seq](0..1, 1/6))
  , axes=none
  , scaling=constrained
):

# In Maple the option orientation is doesn't exist to oriente the text.
# So this trick to suggest you a modification of your answer.

plots:-display(
    plottools:-rotate(pl1, Pi/2)
   , seq(plots:-textplot([0, r, nprintf("%1.2e", r-3819.718635)], align=right), r in [seq](0..6500, 3819.718635/5))
   , seq(
       plots:-textplot(
         [
           7000*cos(t*Pi+Pi/2)
           , 7000*sin(t*Pi+Pi/2)
           , typeset(Pi-t*Pi)
         ]
         , align={left, piecewise(t=1/2, NULL, t > 1/2, below, above)}
       )
       , t in [seq](0..1, 1/6)
     )
);

>

Download suggestion.mw

@delvin Is it this that you want?&n...

Answered: mmcdara 7891

May 21 2023

0 1

@delvin

Is it this that you want?

>

restart;

>	#with(student):

>	U := a[0] + sum(-a[i]tanh(xi[n])^i, i = 1 .. 1) + sum(-b[i]tanh(xi[n])^(-i), i = 1 .. 1):

>	u(xi[n + 1]) := a[0] - a[1](tanh(xi[n]) + tanh(d))/(1 + tanh(xi[n])tanh(d)) - b[1](1 + tanh(xi[n])tanh(d))/(tanh(xi[n]) + tanh(d)):

>	u(xi[n - 1]) := a[0] - a[1](tanh(xi[n]) - tanh(d))/(1 - tanh(xi[n])tanh(d)) - b[1](1 - tanh(xi[n])tanh(d))/(tanh(xi[n]) - tanh(d)):

>	eq := c[1]^2diff(U, xi[n], xi[n])ln(1 + U/alpha) - beta(u(xi[n + 1]) - 2U - u(xi[n - 1])):

>	fin1 := simplify(numer(eq)):

>	fin := simplify(subs(tanh(xi[n]) = Psi, fin1)):

>

# As degree(Psi*ln(Psi), Psi) returns FAIL, the idea is to replace ln(...) by
# something else, let's say H.
#
# Observe that EQ[1] is the identity 0=0

lnPsi := select(has, indets(fin, function), ln)[]:
fin_H := eval(fin, lnPsi=H):
deg   := degree(fin_H, Psi);
print():

for i from 0 to deg do
  collect(coeff(fin_H, Psi, i), H);
  EQ[i] := eval(%, H=lnPsi) = 0;
end do:

print~(EQ):

-2*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))*tanh(d)^2*b[1]*c[1]^2 = 0

(2*tanh(d)^4*b[1]*c[1]^2+2*tanh(d)^2*b[1]*c[1]^2+2*b[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))-2*tanh(d)^2*beta*b[1] = 0

(-2*tanh(d)^4*b[1]*c[1]^2+2*tanh(d)^2*a[1]*c[1]^2-2*tanh(d)^2*b[1]*c[1]^2-2*b[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))+2*tanh(d)^4*beta*b[1]-2*tanh(d)^2*beta*a[1]+2*beta*b[1] = 0

(-2*tanh(d)^4*a[1]*c[1]^2-2*tanh(d)^2*a[1]*c[1]^2+2*tanh(d)^2*b[1]*c[1]^2-2*a[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))+2*tanh(d)^4*beta*a[1]-2*tanh(d)^2*beta*b[1]+2*beta*a[1] = 0

(2*tanh(d)^4*a[1]*c[1]^2+2*tanh(d)^2*a[1]*c[1]^2+2*a[1]*c[1]^2)*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))-2*tanh(d)^2*beta*a[1] = 0

-2*ln(-(Psi^2*a[1]-Psi*alpha-Psi*a[0]+b[1])/(Psi*alpha))*tanh(d)^2*a[1]*c[1]^2 = 0

(1)

>

Download Degree_3_mmcdara.mw

Here is a simpler and corrected version ...

Answered: mmcdara 7891

May 16 2023

0 6

Here is a simpler and corrected version of your code: note the Alert Maplet displays the content of the TextArea component and no longer "test".
But the three times knocking is still here
SimplerButStill3Times.mw

I will keep digging to try and understand where it comes from but I'm afraid it could be some undesired "side effect" between DocumentTools and Maplets packages (?).

This suspicion comes from this observation:
Here is a slight variant of the previous file
SimplerButStill3Times_2.mw
Proceed this way:

write a number in the text area
click on the shortcut to execute the "line" you wrote
click on the line "void"

The Alert Maplet is displayed once
Changing the number and redoing steps 2 and 3 still displays the Alert Maplet once

One way...

Answered: mmcdara 7891

May 14 2023

2 1

It's likely that more concise solutions do exist

repl_U := [seq(z^(-i) = U[q+1-i]/U(z), i = 1..deg_mx)]:
LHS    := collect(eval(lhs(f), repl_U), [seq(U[q+1-i], i = 1..deg_mx)]):
repl_Y := [seq(z^(-i) = Y[q+1-i]/Y(z), i = 1..deg_mx)]:
RHS    := collect(eval(rhs(f), repl_Y), [seq(Y[q+1-i], i = 1..deg_mx)]):
# last operation
LHS = eval(RHS, Y(z)=Y[q+1]);

One_way.mw

In OdeSys you have this ode;diff(Phi(Y),...

Answered: mmcdara 7891

May 13 2023

0 9

In OdeSys you have this ode;

diff(Phi(Y), Y)+.1*Theta(Y)-.1*Phi(Y) = 0;

The boundary condition you set is

Phi(1) = 1+g*D(Phi)(1);

dsolve/numeric doesn't seem to accept a boundary condition of this form for an ode of order one (even if a solution can exist).
See this notional example Download Example.mw

The same problem holds for the ode diff(U(Y), Y) = ...

If you replace your bcs by these one you will get a solution

Cond := 
D(W)(0) = 0, 
W(1) = -delta*(1+1/bet)*D(W)(1), 

D(Theta)(0) = 0, 
Theta(1) = 1+g*D(Theta)(1), 

#U(1) = -delta*(1+1/bet)*D(U)(1), 
U(0)=0,
#Phi(1) = 1+g*D(Phi)(1);
Phi(0)=0

But maybe it's not what you want (see below).

A question for more skilled Maple users: Is it true that dsolve/numeric is not capable to handle a Fourier boundary condition for a first order ode?

ABOUT YOUR FIRST QUESTION : "(I) want extra Boundary conditions U(0)=0 and Phi(0)=0"

A simple idea to achieve this is: instead of solving (for instance)

diff(Phi(Y), Y)+.1*Theta(Y)-.1*Phi(Y) = 0;
Phi(1) = 1+g*D(Phi)(1);

solve instead

diff(diff(Phi(Y), Y)+.1*Theta(Y)-.1*Phi(Y), Y) = 0;
Phi(1) = 1+g*D(Phi)(1);
Phi(0)=0;

Here is an excerpt of the code with a focus on the modifications

#..........

>

OdeSys :=
a1*(1+1/bet)*(diff(W(Y), Y, Y))+a2*Ra*(diff(W(Y), Y))+A-a5*M*W(Y)-S2*W(Y)^2+a2*Gr*Theta(Y)-S*betu*(W(Y)-U(Y)) = 0,
(a4+Rd)*(diff(Theta(Y), Y, Y))+a3*Pr*Ra*(diff(Theta(Y), Y))+Q*Theta(Y)+Pr*alp*S*bett*(Theta(Y)-Phi(Y))+Pr*Ec*((1+1/bet)*a1*(diff(W(Y), Y))^2+a5*M*W(Y)^2+(1+1/bet)*a1*S1*W(Y)^2+S2*W(Y)^3+S*betu*(W(Y)-U(Y))) = 0,
diff(Ra*(diff(U(Y), Y))+betu*(W(Y)-U(Y)), Y) = 0,
diff( Ra*(diff(Phi(Y), Y))+bett*(Theta(Y)-Phi(Y)), Y) = 0

(13)

>	Cond := D(W)(0) = 0, W(1) = -delta(1+1/bet)D(W)(1), D(Theta)(0) = 0, Theta(1) = 1+gD(Theta)(1), U(1) = -delta(1+1/bet)D(U)(1), U(0) = 0, Phi(1) = 1+gD(Phi)(1), Phi(0) = 0;

(14)

>	MVals := [1, 3, 5]

(15)

>	for j to numelems(MVals) do Ans[j] := dsolve(eval([OdeSys, Cond], M = MVals[j]), numeric, output = listprocedure) end do

(16)

>

Full code is here:

Download Worksheet_1D.mw

Solutions...

Answered: mmcdara 7891

May 13 2023

1 0

Two methods are given:

"By hand", which is trivial for problem 1 and needs a little more attention for problem 2 (the 2 body problem).
With Maple package Tolerances.
But I guess that what you are asking for is to prove you can do it yourself? So consider the use of this package as a verification of what you did by hand, or as a tool to solve easily more complex problems you could face.

In addition a more general method based on Taylor expansion is given forthe 2 body problem.

>

restart

A SOLUTION BY HAND

>

local I:

Warning, The imaginary unit, I, has been renamed _I

>	m := (R-delta__R)(I-delta__I)^2(t-delta__t); M := (R+delta__R)(I+delta__I)^2(t+delta__t);

(1)

>	Q__min := eval(m, [R=1.5, delta__R=0.5, I=2.5, delta__I=0.05, t=10, delta__t=0.5]); Q__max := eval(M, [R=1.5, delta__R=0.5, I=2.5, delta__I=0.05, t=10, delta__t=0.5])

(2)

>	Q__mean := (Q__max+Q__min)/2; Q__dev := (Q__max-Q__min)/2; rel_dev := evalf[2](Q__dev/Q__mean)

(3)

>

restart

A READY MADE SOLUTION

>	with(Tolerances):

>

local I:

Warning, The imaginary unit, I, has been renamed _I

>	R := 1.5 &+- 0.5: I := 2.5 &+- 0.05: t := 10 &+- 0.5:

>	Q := RI^2t

(4)

>	Q__mean := NominalValue(Q); Q__dev := ToleranceValue(Q); rel_dev := evalf[2](Q__dev/Q__mean)

(5)

>

restart

THE 2 BODY PROBLEM BY HAND

Remark: the relativae deviation CANNOT be negative!!!

>	F := (y, M, m, R) -> yMm/R^2

(6)

>	F__min := F(y, 7-0.05, 2-0.02, 5+0.05); F__max := F(y, 7+0.05, 2+0.02, 5-0.05);

(7)

>	F__mean := (F__max+F__min)/2; F__dev := (F__max-F__min)/2; rel_dev := evalf[2](F__dev/F__mean)

(8)

AND WITH PACKAGE Tolerances

>

restart

>	with(Tolerances):

>	M := 7 &+- 0.05: m := 2 &+- 0.02: R := 5 &+- 0.05:

>	`F/y` := M*m/R^2

(9)

>	F__mean := NominalValue(`F/y`); F__dev := ToleranceValue(`F/y`); rel_dev := evalf[2](F__dev/F__mean)

(10)

>

restart

THE 2 BODY PROBLEM THROUGH TAYLOR EXPANSION

A more general approach for more general problems.

>	F := yMm/R^2: F__Tayl := mtaylor(F, [M=7, m=2, R=5], 2)

(11)

>	eval(F__Tayl, [M=7+delta__M, m=2+delta__m, R=5+delta__R]): f := unapply(%, [delta__M, delta__m, delta__R])

(12)

>	signs := (signum@diff)~(F, [M, m,R]) assuming M > 0, m > 0, R > 0, y > 0

(13)

>	devs := [0.05, 0.02, 0.05]; F__max := f(op(devs~signs)); F__min := f(op(devs~(-~signs)))

(14)

>	F__mean := (F__max+F__min)/2; F__dev := (F__max-F__min)/2; rel_dev := evalf[2](F__dev/F__mean)

(15)

>

Download Tolerances.mw

To conclude, here is a piece of code aimed to return the upper and lower formal bounds of an arbitrary function F.
Once got you just have to evaluate these bounds with ad hoc nominal valus (mu) and absolute deviations (delta) and compute what you want (the nominal value of F and its total or relative deviation).
You will see that the upper and lower bounds depend upon the order of the Taylor expansion.
MinMaxValues.mw
For example 1 the highest non nul derivative is diff(diff(diff(Q, t), R), I$2) and thus pushing the Taylor expansion further the order 1+1+2=4 is useless.
This is not the case for example 2.

Last but not least, you can use de Units package to get the unit of the absolute deviation and nominal value. Here is an illustration for your first example.
Example1_ByHand_WithUnits.mw

Using series expansion of the res......

Answered: mmcdara 7891

May 12 2023

0 0

Using a series expansion of the rhs term independent of f(t) gets a formal solution in terms wrt t and V.
This could be more interesting than sweeping the V range to explore how the solution behaves:

>

restart

>	U := .21.2894606222-.21(1.158611162/((1+0.4836862955e-1t)(1/(1+0.4836862955e-1*t))^.1986136226)); plot(U, t=0..40)

0.6078673066e-1-.2433083440/((1+0.4836862955e-1*t)*(1/(1+0.4836862955e-1*t))^.1986136226)

>	S := convert(series(0.6078673066e-1-.2433083440/((1+0.4836862955e-1t)(1/(1+0.4836862955e-1t))^.1986136226), t=20), polynom); plot([U, S], t=0..40, color=[blue, red]); ic := f(0) = .21: series_sol := evalf(rhs(dsolve({diff(f(t), t) = (S-.2894606222f(t))/V, ic})));

(1)

>	num_sol := dsolve({diff(f(t), t) = eval((U-.2894606222*f(t))/V, V=2), ic}, numeric): use plots in display( plot(eval(series_sol, V=2), t=0..40, color=red), odeplot(num_sol, [t, f(t)], t=0..40, color=blue) ) end use;

>	plot3d(series_sol, t=0..40, V=0..4, style=surface, orientation=[-53, 62, 14], shading=zhue)

>

Download WithSeriesExpansion.mw

Not at all a simple question....

Answered: mmcdara 7891

May 11 2023

2 3

It all depends what you mean by "strip (=rectangle) packing.
For the a priori simpler case of square packing look here https://en.wikipedia.org/wiki/Square_packing_in_a_square.
Let A lthe side length of the square to pack and a the side length of the packing square: there are no general theoritical results for arbitrary (A, a) values, the number of paxking squares is to be found using ad hoc optimization strategies.

of course if you forbid rotations of the packing square the result is trivial (n=floor((A/a)^2)).

For rectangle packing of a rectangle the things grow in complexity (see "Packing identical rectangles in a rectangle" here https://en.wikipedia.org/wiki/Rectangle_packing#:~:text=Rectangle%20packing%20is%20a%20packing,this%20problem%20have%20been%20studied).

More ressources can be found here https://www.ime.usp.br/~egbirgin/packing/ and in the references herein. In particular http://lobato.it/packing/.

If you forbid rotations of the packing rectangles then, if I'm not mistaken, this piece of code gives you the solution:

inputs
- X and Y are side lengths of the rectangle to pack
- x and are side lengths of the packing strip
outputs
- a list of the form [[n1, n2], n] where:
  - n1 and n2 are the numbers of strips depending on their orientations in the (X, Y) rectangle
  - n is the max of n1 and n2
    This number is a lower bound of the number you are looking for.
    For instance the L-Algorithm for data (49,28,8,3) (http://lobato.it/packing/) gives 57 strips as f(49, 28, 8, 3) = 54.

>

restart:

>	f := proc(X, Y, x, y) local n1, n2: n1 := floor(X/x)floor(Y/y)+floor(max((Y-yfloor(Y/y))/x, 0)): n2 := floor(Y/x)floor(X/y)+floor(max((X-yfloor(X/y))/x, 0)): return [[n1, n2], max(n1, n2)] end proc:

>	f(10, 10, 1, 8)

(1)

>	f(10, 9, 1, 8)

(2)

>	f(10, 7, 3, 3)

(3)

>	f(10, 7, 3, 2)

(4)

>

Download RectanglePacking.mw

In case you would be interested by disk and sphere packing, here is an amazing site http://www.packomania.com/.

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

7891 Reputation

22 Badges

MaplePrimes Activity

These are answers submitted by mmcdara

If I'm not mistaken...

A solution...

Here are two ways:...

Waiting for concurrent proposals......

Obvious solution:...

Didn't I already answered simi...

@Kitonum's answer is likely the simp...

@Rouben Rostamian Hi Rouben, Not...

@delvin Is it this that you want?&n...

Here is a simpler and corrected version ...

One way...

In OdeSys you have this ode;diff(Phi(Y),...

Solutions...

Using series expansion of the res......

Not at all a simple question....

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