mmcdara

3437 Reputation

16 Badges

5 years, 230 days

MaplePrimes Activity


These are answers submitted by mmcdara

An alternative is to use Optimization:-NLPSolve instead of Statistics:-NonLinear.
Both have their own pros and cons.
Here is an example (a few others can be seen in the attached file).
 

restart
X := Vector([1, 2, 3, 4, 5, 6], datatype=float):
Y := Vector([2.2, 3, 4.8, 10.2, 24.5, 75.0], datatype=float):

# pro or con? You must define your objective function

obj := proc(p)
  local f := x -> p[1]+p[2]*x+exp(p[3]*x):
  return add((Y-f~(X))^~2)
end proc:

# admissible intervals for each parameter

lower_bound := `<|>`([6,-5,0]);
upper_bound := `<|>`([7,-4,1]);

Optimization:-NLPSolve(3, obj, [], [lower_bound, upper_bound])


You will find in the attached files several variants of your problem that show why NLPSolve can be a good alternative to NonlinearFit.
But for the problem you submitted  tom's answer seems enough.

reg.mw

@jrive  @acer : This is not engineering notation, 6e-10 should be 600e-12 for instance
difference-between-scientific-engineering-notation-8567943.html

For instance we don't say that 6e-10s equal 0.6 nano-second but 600 pico-second

The procedure EngForm could probably be simplified.

restart

EngForm := proc(x)
  local n, L, p, r, m, e, f:
  if abs(x) < 10^4 and abs(x) > 1 then
    cat(`#mo("`, x, `")`)
  elif abs(x) > 10^4 then
    n := length(op(1, x));
    L := log[10](abs(x));
    p := floor(L);
    r := irem(p, 3);
    m := 10^(L-p+r);
    m := evalf[n](m);
    e := p-r;
    cat(`#mrow(mo("`, m, `"),mo("&#xd7;"),msup(mo("10"),mo("`, e, `")))`)
  elif abs(x) < 1 then
    n := length(op(1, x));
    L := -ceil(log[10](abs(x)));
    r := 3-irem(L, 3);
    p := L+r;
    m := x*10^(p);
    e := -p;
    f := max(r, n);
    m := nprintf(cat("%", r, ".", f-r, "f"), m);
    cat(`#mrow(mo("`, m, `"),mo("&#xd7;"),msup(mo("10"),mo("`, e, `")))`);
  end if;
end proc:
 

EngForm(6e-10);
EngForm(123.4);
EngForm(1.23456e8);
 

`#mrow(mo("600"),mo("&#xd7;"),msup(mo("10"),mo("-12")))`

 

`#mo("123.4")`

 

`#mrow(mo("123.456"),mo("&#xd7;"),msup(mo("10"),mo("6")))`

(1)

Lres := 1/((2*Pi*freq)^2*Cres):

ftest := 10e6:

plot(eval(Lres, freq = ftest), Cres = 0.0 .. 1e-9,

     labels = [Cres, 'Lres'],

     legend = EN(ftest),

     color = red,

     title = 'Inductance*Value*as*a*Function*of*Resonant*Capacitance',

     axis = [gridlines = [default]],
     tickmarks=[[seq(k*10^(-10)=EngForm(k*10^(-10)), k in [seq](2..10, 2))], default],
     size=[700,300]
);

 

 


 

Download EngineerNotation.mw

 

 

 

Same results with a diffferent solving strategy 
 

restart

ode__u__1     := diff(u__1(y),y,y) + 1/2*diff(N(y),y) + 1/2*theta__1(y) = 0;

ode__u__2     := diff(u__2(y),y,y) + theta__2(y) = 0;

ode__N        := diff(N(y),y,y) - 2/3 * (2*N(y) + diff(u__1(y),y)) = 0;

ode__theta__1 := diff(theta__1(y),y,y) = 0;

ode__theta__2 := diff(theta__2(y),y,y) = 0;

diff(diff(u__1(y), y), y)+(1/2)*(diff(N(y), y))+(1/2)*theta__1(y) = 0

 

diff(diff(u__2(y), y), y)+theta__2(y) = 0

 

diff(diff(N(y), y), y)-(4/3)*N(y)-(2/3)*(diff(u__1(y), y)) = 0

 

diff(diff(theta__1(y), y), y) = 0

 

diff(diff(theta__2(y), y), y) = 0

(1)

# start solving ode__theta__1 and ode__theta__2

sol__theta__1 := dsolve({ode__theta__1, theta__1(-1)=1, D(theta__1)(0)=d__theta});
sol__theta__2 := dsolve({ode__theta__2, theta__2(+1)=0, D(theta__2)(0)=d__theta});

theta__1(y) = y*d__theta+d__theta+1

 

theta__2(y) = y*d__theta-d__theta

(2)

# write that theta__1(0) = theta__2(0)

d__theta := solve( eval(eval(theta__1(y) = theta__2(y), [sol__theta__1, sol__theta__2]), y=0) )

-1/2

(3)

# and thus

sol__theta__1 := eval(sol__theta__1);
sol__theta__2 := eval(sol__theta__2);

theta__1(y) = -(1/2)*y+1/2

 

theta__2(y) = -(1/2)*y+1/2

(4)

# which leads to

ode__u__1 := eval(ode__u__1, sol__theta__1);
ode__u__2 := eval(ode__u__2, sol__theta__2);

diff(diff(u__1(y), y), y)+(1/2)*(diff(N(y), y))-(1/4)*y+1/4 = 0

 

diff(diff(u__2(y), y), y)-(1/2)*y+1/2 = 0

(5)

# solve ode__u__2

sol__u__2 := dsolve({ode__u__2, u__2(1)=0, D(u__2)(0)=2*d__u});

u__2(y) = (1/12)*y^3-(1/4)*y^2+2*d__u*y+1/6-2*d__u

(6)

# solve ode__u__1 and ode__N

sols__uN := dsolve({ode__u__1, ode__N, u__1(-1)=0, D(u__1)(0)+N(0)/2=d__u, N(-1)=0, D(N)(0)=0})

{N(y) = (1/12)*exp(-y)*(5+2*exp(-1)+8*d__u)/(exp(1)+exp(-1))-(1/12)*exp(y)*(-5-8*d__u+2*exp(1))/(exp(1)+exp(-1))-1/6-(1/12)*y^2+(1/6)*y-(2/3)*d__u, u__1(y) = (1/24)*exp(-y)*(5+2*exp(-1)+8*d__u)/(exp(1)+exp(-1))+(1/24)*exp(y)*(-5-8*d__u+2*exp(1))/(exp(1)+exp(-1))+(1/18)*y^3-(1/6)*y^2+(1/12+(4/3)*d__u)*y-(1/72)*(12*exp(1)*exp(-1)-72*exp(1)*d__u-120*exp(-1)*d__u-7*exp(1)-37*exp(-1))/(exp(1)+exp(-1))}

(7)

# write the conditions u__1(0)=u__2(0)

d__u := solve( eval(eval(u__1(y)=u__2(y), sols__uN union {sol__u__2}), y=0) )

(1/24)*(12*exp(1)*exp(-1)-exp(1)-31*exp(-1))/(9*exp(1)+11*exp(-1))

(8)

sols := eval({sols__uN[], sol__u__2,  sol__theta__1, sol__theta__2}):
print~(simplify(sols)):

N(y) = -(1/36)*(27*exp(2+y)*y^2+54*exp(2+2*y)-54*exp(2+y)*y+53*exp(2+y)-134*exp(2*y+1)+33*y^2*exp(y)+12*exp(1+y)-66*y*exp(y)-134*exp(1)+35*exp(y)-66)*exp(-y)/(9*exp(2)+11)

 

u__1(y) = (1/36)*(18*exp(2+y)*y^3-54*exp(2+y)*y^2+27*exp(2+2*y)+25*exp(2+y)*y+22*y^3*exp(y)+30*exp(2+y)-67*exp(2*y+1)+24*exp(1+y)*y-66*y^2*exp(y)-36*exp(1+y)-29*y*exp(y)+67*exp(1)+126*exp(y)+33)*exp(-y)/(9*exp(2)+11)

 

u__2(y) = (1/12)*(9*exp(2)*y^3-27*exp(2)*y^2-exp(2)*y+11*y^3+19*exp(2)+12*exp(1)*y-33*y^2-12*exp(1)-31*y+53)/(9*exp(2)+11)

 

theta__1(y) = -(1/2)*y+1/2

 

theta__2(y) = -(1/2)*y+1/2

(9)

plots:-display(
  plot(eval(u__1(y), sols), y=-1..0, color=red),
  plot(eval(u__2(y), sols), y=0..1, color=blue)
);

plots:-display(
  plot(eval(u__1(y), sols), y=-0.01..0, color=red),
  plot(eval(u__2(y), sols), y=0..0.01, color=blue)
)

 

 

plot(eval(N(y), sols), y=-1..0)

 

print~(simplify(sols)):

N(y) = -(1/36)*(27*exp(2+y)*y^2+54*exp(2+2*y)-54*exp(2+y)*y+53*exp(2+y)-134*exp(2*y+1)+33*y^2*exp(y)+12*exp(1+y)-66*y*exp(y)-134*exp(1)+35*exp(y)-66)*exp(-y)/(9*exp(2)+11)

 

u__1(y) = (1/36)*(18*exp(2+y)*y^3-54*exp(2+y)*y^2+27*exp(2+2*y)+25*exp(2+y)*y+22*y^3*exp(y)+30*exp(2+y)-67*exp(2*y+1)+24*exp(1+y)*y-66*y^2*exp(y)-36*exp(1+y)-29*y*exp(y)+67*exp(1)+126*exp(y)+33)*exp(-y)/(9*exp(2)+11)

 

u__2(y) = (1/12)*(9*exp(2)*y^3-27*exp(2)*y^2-exp(2)*y+11*y^3+19*exp(2)+12*exp(1)*y-33*y^2-12*exp(1)-31*y+53)/(9*exp(2)+11)

 

theta__1(y) = -(1/2)*y+1/2

 

theta__2(y) = -(1/2)*y+1/2

(10)

 


 

Download 3PBVP.mw

 

A way.
Note that a Haar function should depend on 2 parameters (j and k in the attached file).

 

restart:

# Normalized Haar functions
# In your case replace 2^(j/2) by 1 and -2^(j/2) by -1
#

Haar := (j, k) -> t -> `if`(
                          j = 0,
                          piecewise(0 <= t and t <= 1, 1, 0),
                          piecewise(
                            (k-1)/2^j   <= t and t <= (k-1/2)/2^j,  2^(j/2),
                            (k-1/2)/2^j <= t and t <=  k/2^j     , -2^(j/2),
                                                                    0
                          )
                       );

proc (j, k) options operator, arrow; proc (t) options operator, arrow; `if`(j = 0, piecewise(0 <= t and t <= 1, 1, 0), piecewise((k-1)/2^j <= t and t <= (k-1/2)/2^j, 2^((1/2)*j), (k-1/2)/2^j <= t and t <= k/2^j, -2^((1/2)*j), 0)) end proc end proc

(1)

plot(Haar(0, 0)(t), t=0..1)

 

# I use a vertical shift (d*k) for a better visualization

n := 2;
d := 2*(2^(n/2)):
plot(
  [seq(d*k + Haar(n, k)(t), k=1..2^n)],
  t=0..1,
  legend=[seq([n, k], k=1..2^n)],
  color=[seq(ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]), k=1..2^n)],
  axis[2]=[tickmarks=[seq(d*k=0, k=1..2^n)]],
  gridlines=true
)

2

 

 

n := 3;
d := 2*(2^(n/2)):
plot(
  [seq(d*k + Haar(n, k)(t), k=1..2^n)],
  t=0..1,
  legend=[seq([n, k], k=1..2^n)],
  color=[seq(ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]), k=1..2^n)],
  axis[2]=[tickmarks=[seq(d*k=0, k=1..2^n)]],
  gridlines=true
)

3

 

 

 


 

Download Haar.mw



this could be done by hand without any difficulty and using Maple for this is seems unnecessary
Nevertheless, look to the answers in the attached file.

PS 1:  for question 2 the values of p are complex
PS 2:  maybe the "trouble you have" in solving this comes from the fact that you omited the multiplicative symbol almost everywhere in the expressions Question 2 contains???
 

A sledgehammer to kill a fly   Download Q1_Q2.mw

If your homework requires more handy work look here  Download Q1_elem.mw


 

restart;

kernelopts(version);

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

expr := -.27059805007*sin(.120000+epsilon)

        +.27059805007000*sqrt(1.-cos(.12+epsilon)^2);

-.27059805007*sin(.120000+epsilon)+.27059805007000*(1.-cos(.12+epsilon)^2)^(1/2)

(2)

angles := op~([indets(expr, function)[]])

[.12+epsilon, .120000+epsilon]

(3)

Expr := eval(expr, angles =~ [z$2]):
simplify(Expr, trig) assuming z in RealRange(0, Pi)

0.

(4)

 


 

Download Maple2015.mw


you would like to r

  • rewrite a differential equation, let's say
  • a*diff(f(x), x)+b*f(x)=c

in a symbolic form 

  • (a*D + b*I)(f) = c

where I is the identity operator and D the differential one.

  • Next to define the formal inverse L of (a*D + b*I)
    L = (a*D + b*I)-1
    
  • And use L to simplify another ODE which contains f(x) by replacing each of its occurrences by L(c)


This is possible only in very specific situations (see at the end of the attached file)
This file unfolds the previous "program" step by step and shows its limits when placed in a general framework.

Operator_manipulations.mw

Thus, the answer to your question "can it be done by factoring differential operators and NOT directly using dsolve?" is unfortunatele NO (or it is far beyond my capabilities)


Identifying a list of integers to a random variable may be  necessary or not (for instance tandomly chose the rows or columns of a matrix, as it happens in some numerical algorithms).
As your list is not a list of contiguous or evenly space elements, maybe the statistical approach can be useful?

Here are a few alternatives based on an identification to a random variable (the last alternative doesn't).
So make your choice.
 

restart;

 

Alternative 1a

 

L := [ seq(-720..-30, 15), seq(360..720, 15) ]:

with(Statistics):

# This step is necessary to get integer samples.
# The default UseHardwareFloats := true forces Sample
# to return Hfloats makink round unefficient.

UseHardwareFloats := false:

# Define a discrete uniform random variable whose realizations are in L

X := RandomVariable(EmpiricalDistribution(L)):

# Draw one single sample from X
round(Sample(X, 1)[1]);

# or several samples
round~(Sample(X, 10));

-660

 

Vector[row](%id = 18446744078296578278)

(1)

 

Alternative 1b

 

restart;

L := [ seq(-720..-30, 15), seq(360..720, 15) ]:

with(Statistics):

# Shuffle makes a random permutation of L, one just need to take the first element of
# the shuffled list.
#
# This alternative is very closed to "Alternative 3" above

Shuffle(L)[1]

-225

(2)

# Works with non numeric lists too:


L := ["flower", "cat", "mouse", "sun"]:
Shuffle(L)[1]

"mouse"

(3)

 

Alternative 2

 

restart;

L := [ seq(-720..-30, 15), seq(360..720, 15) ]:

with(Student[Statistics]):

X := EmpiricalRandomVariable(L):

# Much simpler than Statistics for Sample really samples X,
# not some Hloat representation of it.
#
# This is my favorite alternative

Sample(X, 1);

# or

Sample(X, 10)

Vector[row](1, {(1) = 525})

 

Vector[row](%id = 18446744078300839934)

(4)

 

Alternative 3

 

restart

L := [ seq(-720..-30, 15), seq(360..720, 15) ]:

with(combinat):

# The first term of random permutation of L is the desired sample.
#
# Watchout: this alternative is not efficient if the list is large

combinat:-randperm(L)[1];

# or

seq(combinat:-randperm(L)[1], k=1..10);

-420

 

510, 720, -450, 720, -405, -255, -285, 525, -210, 690

(5)

# A hammer to kill a fly?
# But it can be usefull for non numeric lists, see "Alternative 1b" above
# (a thing tomleslie's proposal can do too):

L := ["flower", "cat", "mouse", "sun"]:
combinat:-randperm(L)[1];
# or

seq(combinat:-randperm(L)[1], k=1..5);

"flower"

 

"mouse", "flower", "cat", "cat", "mouse"

(6)

 


 

Download Alternatives.mw

 


One thing: you can set a=1 without any loss of generality (same if you set d=1).

With this in mind it becomes extremely easy to construct by hand (almost by hand) particular solutions.
In the attached file one respects the constraints that the coefficients are in [1, 10], another don't.

galimatias.mw

Of course it's very far from Kitonum's exhaustive search

Look here view.aspx (paragraph "The External C Compiler").

Contrary to Windows platforms
"On all UNIX and Mac OS X systems, the GNU C compiler is not distributed with Maple, and must be installed separately, if it is not already present on the system " 

This could be the cause of your problem

 

One way which introduces only a slight modification to your code is to replace the line

eq1 := Qbar = InvT . Q . R . T . InvR

by the red line below

assign([entries(Qbar, nolist)] =~ [entries(InvT . Q . R . T . InvR, nolist)]);

plot(Qb11, p = 0 .. 9): # fails for Qb11 depends on 3 variables:

Warning, expecting only range variable p in expression (.2339357430e12*cos(p)^2+.2319277108e11*s^2)*cos(p)^2+(.2319277108e11*cos(p)^2+s^2*Q22)*s^2+.2868e11*s^2*cos(p)^2 to be plotted but found names [Q22, s]

# check this
indets(Qb11, name);

                          {Q22, p, s}

# do this, for instance, to evaluate Qb11 for siome values of Q22 and s
plot(eval(Qb11, [Q22=1, s=1]), p=0..9)

 

with recent Maple versions

numer(a) %/ denom(a);


For older versions (2015 for instance)
 

with(InertForm):
n := convert(numer(a), string):
d := convert(denom(a), string):
Parse( cat("(", n, ")/", d) ):
Display(%);

 


Hints for Problem 7
An equivalence relation E has 3 properties R, S and T (which are?).

  • one of them (R) is obvious
  • one (S) comes from the commutativity of the addition
  • the last (T) involes 3 couples C, C' and C"
    • write C E C' and C' E C""
    • try to combine the two equalities to obtain C E C" (needs a simplification and uses properties of the addition)
       
  • In fact, let F a relation which has the same properties of the addition, then  (a, b) E (c, d) defined by  (a+d) F (b+c) is an equivalence relation.

 

For the three other problems look at some hints here
pbs.pdf


It's clear that D is at the intersection of the bisector of angle A and the line which passes through B and C.
Thus its coordinates and a parametric representation of AD

coords__D := ([-6, 1,-1] +~ [2,-3, 7]) /~ 2;
line__AD  := [x, y, z] =~ [-2, 3, 5]*~s +~ coords__D*~(1-s);

 

But as I'm fond of the geom3d package, here is a step by step wauy yo use it (uselesss for this problem but it could help you for more complex ones). 

restart:
with(geom3d):

point(A, -2, 3, 5):
point(B, -6, 1,-1):
point(C,  2,-3, 7):

line(AB, [A, B]):
line(AC, [A, C]):

# parametric equations of AB and AC

eq__AB := [x, y, z] =~ Equation(AB, s):
eq__AC := [x, y, z] =~ Equation(AC, s):

# AD is the bisector of angle BAC:
# parametric equation of AD

eq__AD := (eq__AB +~ eq__AC) /~ 2;
               [x = -2, y = 3 - 4 s, z = 5 - 2 s]

# implicit equation of AD 
map(u -> if depends(u, s) then solve(u, s) else rhs(u) end if, eq__AD):
eq__AD__imp := add(%);
                         5   1     1  
                         - - - y - - z
                         4   4     2  

# coordinates of D
local D:
line(BC, [B, C]):
line(AD, rhs~(eq__AD), s):
intersection(D, BC, AD):
detail(D)
          Geom3dDetail(["name of the object", D], 
            ["form of the object", point3d], 
            ["coordinates of the point", [-2, -1, 3]])

# visualize
plots:-display(
  PLOT3D(
    POINTS([coordinates(A)], SYMBOL(_SOLIDCIRCLE, 30), COLOR(RGB, 1, 0, 0)),
    POINTS([coordinates(B)], SYMBOL(_SOLIDCIRCLE, 30), COLOR(RGB, 1, 0, 0)),
    POINTS([coordinates(C)], SYMBOL(_SOLIDCIRCLE, 30), COLOR(RGB, 1, 0, 0)),
    POINTS([coordinates(D)], SYMBOL(_SOLIDCIRCLE, 30), COLOR(RGB, 0, 0, 1)),
    TEXT(1.1*~coordinates(A), "A", FONT(Helvetica, bold, 14)),
    TEXT(1.1*~coordinates(B), "B", FONT(Helvetica, bold, 14)),
    TEXT(1.1*~coordinates(C), "C", FONT(Helvetica, bold, 14)),
    TEXT(1.2*~coordinates(D), "D", FONT(Helvetica, bold, 14), COLOR(RGB, 0, 0, 1)),
    POLYGONS([coordinates(A), coordinates(B), coordinates(C)], COLOR(RGB, 0.9$3))
  ),
  plot3d([rhs~(eq__AD)[]], s=0..1, color=blue, thickness=5)
)



with_geom3d.mw

I'm not going to risk a stiff neck to read this but it sounds like a fixed point method.
Here is an example, adapt it to your own function f


 

restart:

with(plots):

# example

f := x -> surd(x, 3):
g := x -> surd(x, 4) / 2:

pf := plot(f(x), x=0..1.3, color=blue , thickness=2, gridlines=true):
pg := plot(g(x), x=0..1.3, color=green, thickness=2, gridlines=true):
pl := plot( x  , x=0..1.3, color=black):

start := eval(x, solve({g(x)=x, x > 0}));

(1/4)*2^(2/3)

(1)

tu := `#mo("&#x25b2;")`;
tr := `#mo("&#x25ba;")`;

`#mo("&#x25b2;")`

 

`#mo("&#x25ba;")`

(2)

iter := proc(p, eps)
  local q := copy(p):
  local L, r:
  L := q, [q[1], f(q[1])];
  r := [f(q[1]), f(q[1])];
  L := L, r;
  while add(q-~r)^~2 > eps^2 do
    q := r:
    L := L, [q[1], f(q[1])];
    r := [f(q[1]), f(q[1])];
    L := L, r;
  end do:
  return L
end proc:

path := [iter([evalf(start), 0], 0.05)]:
display(
  pf, pg, pl,
  plot(path, linestyle=3),
  seq(textplot([op((path[k]+~path[k+1])/~2), tr]), k in [seq(2..numelems(path), 2)]),
  textplot([op(([start$2]+~path[2])/~2), tu]),
  seq(textplot([op((path[k]+~path[k+1])/~2), tu]), k in [seq(3..numelems(path)-1, 2)]),
  size=[800, 800]
)

 

 


 

Download fixed_point.mw

Look to the first line in the attached fie, this could also interest you

 

4 5 6 7 8 9 10 Last Page 6 of 28