mmcdara

3713 Reputation

17 Badges

5 years, 347 days

MaplePrimes Activity


These are answers submitted by mmcdara


In fact the true question is "What form would you like this sensitivity report to take?"

Here are a few solutions (I used a larger step to provide "small" displays.
The result matrix can be saved using ExportMatrix for instance.

SR.mw

Unless I made a mistake here is a proof
 

restart:

g := exp(x)+exp(y)+exp(z)+2*exp(-x-y-z)-4*2^(1/4);

exp(x)+exp(y)+exp(z)+2*exp(-x-y-z)-4*2^(1/4)

(1)

# Write
#   u=exp(x)
#   v=exp(y)
#   w=exp(z)

h := u+v+w+2/(u*v*w)-4*2^(1/4)

u+v+w+2/(u*v*w)-4*2^(1/4)

(2)

# Observe that:
#   u+v+w is 3 times the Arithmetic Mean (AM) of (u, v, w)
#   u*v*w equals the Geometric Mean (GM) of (u, v, w) raised to the cube
# Thus

h := 3*AM+2/GM^3-4*2^(1/4)

3*AM+2/GM^3-4*2^(1/4)

(3)

# Let's focus on the term 3*AM+2/GM^3
#
# From the AM–GM inequality
# https://en.wikipedia.org/wiki/Inequality_of_arithmetic_and_geometric_means
# one has GM <= AM
# Thus:
#   GM^3 <= AM^3
#   1/GM^3 >= 1/AM^3
#   3*AM+2/GM^3 >= 3*AM+2/AM^3
#
# Solve the equation 3*AM+2/AM^3 > 4*2^(1/4) given that AM > 0

solve(3*AM+2/AM^3 > 4*2^(1/4)) assuming AM > 0

RealRange(Open(0), Open(2^(1/4))), RealRange(Open(2^(1/4)), infinity)

(4)

# Thus 3*AM+2/AM^3 is
#    positive provided AM <>2^1/4,
#    null if AM = 2^1/4,
# and, because AM is necessarilly > 0 (arithemic mean of exponentials),
# this proves 3*AM+2/AM^3 > 4*2^(1/4) is true but at AM = 2^1/4.
#
# Given that 3*AM+2/GM^3 >= 3*AM+2/AM^3 you have your proof.

 

Download AM_GM.mw

An even simpler way based on th same AM-GM inequality is
 

# Other way: AM = GM iif u=v=w.
# Then it's enough to demonstrate that 3*u+2/u^3 cannot be less than 4*2^(1/4)
# for any strictly positive u:

infolevel[solve]:=1:
solve({3*u+2/u^3 < 4*2^(1/4), u > 0});
solve: Warning: no solutions found

 

Possibly apart the obvious loops at edges x0 and x10(*) I can't see any cycle (more precisely circuit as your graph is directed) on your graph !
A cycle circuit in a directed graph G is a non-empty path in G with its first vertex equal to its last vertex.

(*) some will consider that a loop is not circuit; ths choice to do this or not generally depends upon the context

If I correctly understood your problem you will find several ways to do this in the attached file.
 

restart:

One way is to define a function with 4 indeterminates

f   := (a, b, c, x )-> a*x^2+b*x+c;
sol := solve(f(a, b, c, x)=0, x);

# ... the following is strictly identical to what id done below.

 

proc (a, b, c, x) options operator, arrow; a*x^2+b*x+c end proc

 

(1/2)*(-b+(-4*a*c+b^2)^(1/2))/a, -(1/2)*(b+(-4*a*c+b^2)^(1/2))/a

(1)

Another way is to define a parameterized function of x this way

f := (a, b, c) -> x -> a*x^2+b*x+c;

proc (a, b, c) options operator, arrow; proc (x) options operator, arrow; a*x^2+b*x+c end proc end proc

(2)

sol   := solve(f(a, b, c)(x)=0, x):
optx  := unapply( sol[1], (a, b, c));
opttx := unapply( sol[2], (a, b, c));
 

proc (a, b, c) options operator, arrow; (1/2)*(-b+(-4*a*c+b^2)^(1/2))/a end proc

 

proc (a, b, c) options operator, arrow; -(1/2)*(b+(-4*a*c+b^2)^(1/2))/a end proc

(3)

# to evaluate optx at b=2 and c=4 just do

optx(a, 2, 3);
 

(1/2)*(-2+(-12*a+4)^(1/2))/a

(4)

# be careful, this value can be complex

eval(optx(a, 2, 3), a=-1);
evalf(%);

eval(optx(a, 2, 3), a=1);
evalf(%);

1-(1/2)*16^(1/2)

 

-1.000000000

 

-1+(1/2)*(-8)^(1/2)

 

-1.+1.414213562*I

(5)

# A plot without caution (Maple 2015)
# Note that the plot is right for a=-1, but wrong for a=+1

plot(optx(a, 2, 3), a=-2..2, gridlines=true);

 

# A more convenient way

plot([Re(optx(a, 2, 3)), Im(optx(a, 2, 3))], a=-2..2, color=[blue, red], legend=['Re(optx)', 'Im(optx)'], gridlines=true);

 

A third way could be to define this parameterized function as below, but it seems less natural to me
(my proper opinion)

f   := x -> (a, b, c) -> a*x^2+b*x+c;
sol := solve(f(x)(a, b, c)=0, x);

proc (x) options operator, arrow; proc (a, b, c) options operator, arrow; a*x^2+b*x+c end proc end proc

 

(1/2)*(-b+(-4*a*c+b^2)^(1/2))/a, -(1/2)*(b+(-4*a*c+b^2)^(1/2))/a

(6)

 


 

Download Several_ways.mw

This just says that the solution of F(T, A)=0 with respect to T expresses itself in terms of the LambertW function.

Your solution cannot be expressed in another way given the definition of the LambertW function.
(In fact this latter has been introduced to designate the series solution of the the equation you want to solve; see
here Lambert_W_function (History) for more details)
 

restart

df := (-7.233333422*10^10*exp(.45*T)+1.617777797*10^11)/T-(-1.607407427*10^11*exp(.45*T)+1.617777797*10^11*T+1.745679033*10^11+A)/T^2

(-0.7233333422e11*exp(.45*T)+0.1617777797e12)/T-(-0.1607407427e12*exp(.45*T)+0.1617777797e12*T+0.1745679033e12+A)/T^2

(1)

solve(df, T);

2.222222222*LambertW(-.3995249844-0.2288650873e-11*A)+2.222222222

(2)

# What does LambertW mean?

#help(LambertW)

# Example

exple := (a*x+b)*exp(k*x)=c;;
solve(exple, x);

(a*x+b)*exp(k*x) = c

 

(a*LambertW(c*k*exp(b*k/a)/a)-b*k)/(a*k)

(3)

# Assuming T <> 0, df=0 if numer(df)=0

ndf := numer(df);

-0.7233333422e11*T*exp(.45*T)-0.1745679033e12+0.1607407427e12*exp(.45*T)-A

(4)

# Put this into a more convenient form

expr := subs(exp(.45*T)=freeze(exp(.45*T)), ndf):
e    := indets(%, name)[-1]:
conv := thaw(collect(expr, e))

(-0.7233333422e11*T+0.1607407427e12)*exp(.45*T)-0.1745679033e12-1.*A

(5)

# "conv" has the same form of "exple" with
#  a = -7.233333422*10^10
#  b = 1.607407427*10^11
#  c = 1.745679033*10^11+1.*A
#  k = .45
#
# thus the solution you get.

 


 

Download LambertW.mw

 


The term probability function for a discrete random variable (RV)  should be conveniently replaced by the term probability mass function (PMF).
The PMF of a discrete RV defined over a countable set S of integers (for instance) is null evertwhere but at the points in S.

I consider IMO that the plots below are more correct.

 

restart;

with(Statistics):

X := RandomVariable(Binomial(45, 0.9)):
P := seq(ProbabilityFunction(X, i), i=0..45):

plots:-display(
  seq(
    plot([ [i, 0], [i, P[i+1]] ], thickness=5),
    i=1..45
  )
);

 

reduced_P := map(n -> if P[n+1] > max([P])/1000 then
                    plot([ [n, 0], [n, P[n+1]] ], thickness=7)
                    end if,
                    [$0..45]
                ):
plots:-display(reduced_P)

 

 


Download Probability_plot.mw


PS,
Let X a discrete RV.
Still IMO I don't agree with Maple returning a non zero value for ProbabilityFunction(X, t) when t doesn't belong to set S.
If this results was good, then the CDF of X wouldn't be a constant in the interval (s[i], s[i+1]) where s[i] and s[i+1] are two consecutive points in S (nor even a discontinuous function).
Use this to convince yourself that Maple is not consistent with the definitions it uses.

# Integrate the probability function between 0 and s
# One might expect this is the CDF of X

F := int(ProbabilityFunction(X, t), t=0..s);

# Compare the plots of F and ths CDF of X

plots:-display(
  plot(F, s=30..46, color=red),
  plot(CDF(X, s), s=30..46, color=blue)
)

 

 This generates a random polynomial whose coefficients are randomly choosen in the list L.

L := [1, 1, 2, 5, 6]:
randpoly(z, coeffs = proc() combinat:-randperm(L)[1] end proc);

 


EDITED 

I hadn't correctly understood your question, sorry.
There is no difficilty at all to do this with Explore:

 

f := proc(a)
plots:-display(
  plot(5, theta = 0 .. 2*Pi, coords = polar, color = black, thickness = 3, filled = [color = yellow, transparency = .9]),
  plot(
    5*sin(a*theta), 
    theta = 0 .. 2*Pi, 
    coords = polar, 
    color = magenta, 
    thickness = 3, 
    filled = [color = green, transparency = .5]
  ),
  plot(
    5*cos(a*theta), 
    theta = 0 .. 2*Pi, 
    coords = polar, 
    color = red, 
    thickness = 3, 
    filled = [color = blue, transparency = .5]
  )
)
end proc:


Explore(f(a), parameters = [a = 2 .. 6])


Download PétalesGraphes_Maple1015.mw

Point 1
What quantity has dimension 1/time^2 : no one, but it shoulf not be a problem
(of course an acceleration divided by a mass is homogenous tio the inverse of  time squared but I don't think its the answer you wnanted).

Point 2
By integration both sides of DIV from 0 and 3 one can avoid computing an integral to find the length of the wire:

a:=0.36: 
DIV := diff(y(x), x, x) = a*sqrt(1 + diff(y(x), x)^2):
RV  := y(0)=2.6, y(3)=2.1:
sol := rhs(dsolve({DIV, RV}, y(x))):
L   := evalf( ( eval(diff(sol, x), x=3) - eval(diff(sol, x), x=0) ) / a):
evalf(value(%));

                          3.187401749


 

restart:

DIV := diff(y(x), x, x) = a*sqrt(1 + diff(y(x), x)^2):
RV  := y(0)=2.6, y(3)=2.1:
sol := rhs(dsolve({DIV, RV}, y(x))):

L    := evalf( ( eval(diff(sol, x), x=3) - eval(diff(sol, x), x=0) ) / a):
vals := evalf([allvalues(L)]):  # 2 opposite solutions
sag  := unapply(vals, a):

plot(
  [abs(sag(a)[1]), abs(sag(a)[2])],
  a=0.2..0.8,
  color=[red, blue],
  thickness=[1, 11],
  transparency=[0, 0.8],
  gridlines=true,
  labels=['a', 'length']
)

 

# which value of a gives a length of 3.6?

rel := a = select(is, fsolve~(vals=~3.6), positive)[]

a = .6900379205

(1)

# Value of the Sag defined as the maximum vertical distance
# between the shape of the wire and the straitght line
# which passes by through the points (0, 2.6) and (3, 3.1):

sola  := eval(sol, rel):
vdist := sola - (2.6+(2.1-2.6)*x/3):

Sag := -minimize(vdist, x=0..3, location);

.8569225658, -{[{x = 0.}, -0.80202e-8], [{x = 1.462158429}, -.8569225658], [{x = 3.}, -0.8651e-8]}

(2)

# thus the solution

Sag__xy := op(2, [Sag][2])[2]

[{x = 1.462158429}, -.8569225658]

(3)

 


Download catenary_4.mw

Alternative to minimize:

x = select(is, [solve(diff(vdist, x)=0, x)], positive)[]:
allvalues([eval(vdist, %)]):
y = select(is, op~([%]), real)[]:

Sag__xy= [%%%, %];
         Sag__xy = [x = 1.462158429, y = -0.8569225635]

 

As Ì said in my answer to tyour previous question, this can be done in a simpler way

restart;
# if you use integer values, I advice you to scale the equations 
# by using a scaling factor C (this will help fsolve finding the solution).
# If you use floats instead this is not necessary

C := 100:

DIV := diff(y(x), x, x) = a*sqrt(1 + diff(y(x), x)^2):
RV  := y(0)=50/C, y(d)=70/C:

sol := rhs(dsolve({DIV, RV}, y(x))):
L   := Int(sqrt(1+diff(sol, x)^2), x=0..d) assuming x>0, d>0:

# abscissa of the contact to the ground

c := solve(diff(sol, x)=0, x):
# equations to find a and d:

eq1 := eval(sol, x=c)=0:
eq2 := value(L)=140/C:
ad  := fsolve({eq1, eq2}, {a, d})

              {a = 9.213375344, d = 0.5541723235}

# unscaled a and d

AD := [a, d]=~eval([a/C, d*C], ad)
              [a = 0.09213375344, d = 55.41723235]

 

@brian bovril 

The boundary conditions you used are not correct; moreover @vv  solution wich consists in solving two half wires seems a little bit artificial and complicated (IMO) .
If the wire has no transverse stiffness (whch is the basic assumptions used to derive the catenary equation), 
then the boundary conditions are 

y(0) = h0, y(d)=h1

where d is the distance between the two poles.

Here is a partial numeric solution for two poles of heights 1 and 2 and a cable of length 4 (just an example to adapt to the true values [see foot note])
 

restart;

# A simpler problem:
#
# Find the distance between two poles of height 1 and 2 such
# that the wire has length 4.

DIV := diff(y(x), x, x) = sqrt(1 + diff(y(x), x)^2):
RV  := y(0)=1, y(d)=2:
sol := rhs(dsolve({DIV, RV}, y(x))):
L   := Int(sqrt(1+diff(sol, x)^2), x=0..d) assuming x>0, d>0:
gap := solve(value(L)=4, d);

plot(eval(sol, d=gap), x=0..gap)

ln(17/2+(1/2)*285^(1/2))

 

 

# As this problem has a single solution we need an extra parameter
# to find this more complicated problem:
#
# Find the distance between two poles of height 1 and 2 such
# that the wire has length 4 AND the wire tangents the ground.

DIV := diff(y(x), x, x) = a*sqrt(1 + diff(y(x), x)^2):
RV  := y(0)=1, y(d)=2:
sol := rhs(dsolve({DIV, RV}, y(x)));
L   := Int(sqrt(1+diff(sol, x)^2), x=0..d) assuming x>0, d>0;

cosh(a*x+ln(RootOf(((exp(a*d))^2-exp(a*d))*_Z^2-2*a*exp(a*d)*_Z-exp(a*d)+1)))/a+(1/2)*(2*a*(exp(a*d))^2*RootOf(((exp(a*d))^2-exp(a*d))*_Z^2-2*a*exp(a*d)*_Z-exp(a*d)+1)-4*a*exp(a*d)*RootOf(((exp(a*d))^2-exp(a*d))*_Z^2-2*a*exp(a*d)*_Z-exp(a*d)+1)-(exp(a*d))^2+1)/(exp(a*d)*(exp(a*d)-1)*a*RootOf(((exp(a*d))^2-exp(a*d))*_Z^2-2*a*exp(a*d)*_Z-exp(a*d)+1))

 

Int((1+sinh(a*x+ln(RootOf(((exp(a*d))^2-exp(a*d))*_Z^2-2*a*exp(a*d)*_Z-exp(a*d)+1)))^2)^(1/2), x = 0 .. d)

(1)

# abscissa of the contact to the ground

c := solve(diff(sol, x)=0, x)

-ln(RootOf(((exp(a*d))^2-exp(a*d))*_Z^2-2*a*exp(a*d)*_Z-exp(a*d)+1))/a

(2)

# equations to find a and d:

eq1 := eval(sol, x=c)=0:
eq2 := value(L)=4:
ad  := fsolve({eq1, eq2}, {a, d})

{a = 1.691767707, d = 2.248911615}

(3)

plot(eval(sol, ad), x=0..eval(d, ad))

 

 

 

Download catenary_2.mw


foot note : with h1=50, h2=70 and L=160 fsolve can't find a solution (I killed it after one minute).
I suggest you to make the problem dimensionless to make it easier to solve.
For instance, let C > 0, and denote

(a, d) = solution of the problem y(0)=h0, y(d)=h1, length=L

then the solution of the problem

y(0)=C*h0, y(d)=C*h1, length=L*C

is 

(a/C, d*C)



Last remark: In a practical problem, a would be an input of the problem (a depends on the properties of the wire: lineic mass and Young modulus) and the length L of the cable would be an output, not an input.

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

5 6 7 8 9 10 11 Last Page 7 of 30