acer

30152 Reputation

29 Badges

18 years, 232 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

For some values of parameter Cl the dsolve command may here return a result containing an inert integral, rather than a fully explicit answer.

In that case you could apply the value command. (This aspect of dsolve is mentioned in a bullet point a section of the ?dsolve,details Help page.)

Eg,

     plot(value(sol(t)), t = 0 .. 48)

dsH_ex.mw

note: There are other ways to deal with it (changing the Int integral to use some specific numeric methods, or doing a numeric dsolve, etc.) Let us know if you'd like to see such.

@jrive First I'll show a form similar to that in nm's Answer, and discuss more on how that is not equivalent to the original even if all parameters are real. (See also nm's pertinent note, "sqrt(a)/b is same as sqrt(a/b^2) assuming b>0", by which he qualifies his answer.)

restart;

sol := (-v + sqrt(-4*a^2*R^2 + v^2))/(2*a*omega*L);

(1/2)*(-v+(-4*R^2*a^2+v^2)^(1/2))/(a*omega*L)

alt := evalindets(expand(sol),`*`,
                  u->`if`(membertype(anything^(1/2),u),sqrt(u^2),u));

-(1/2)*v/(a*omega*L)+(1/2)*((-4*R^2*a^2+v^2)/(a^2*omega^2*L^2))^(1/2)


Example for which they're not equal.

eval([sol,alt], [R=1,v=1,a=1,L=1,omega=-1]);

[1/2-(1/2)*(-3)^(1/2), 1/2+(1/2)*(-3)^(1/2)]

combine(simplify( sol-alt )) assuming real, a*omega*L>0;

0

combine(simplify( sol-alt )) assuming real, a*omega*L<0;

(-4*R^2*a^2+v^2)^(1/2)/(a*omega*L)

Download rad_comb_ex.mw

Regarding your comment about bringing a term into a radical, consider this simpler example,

ee := sqrt(x^2)/x;

(x^2)^(1/2)/x

simplify(ee) assuming x>0;

1

simplify(ee) assuming x<0;

-1


If we absorb that x^1 into the radical (by squaring it) and make no other changes then we've mistakenly lost the multiplicative sign of x. If x is positive then it's ok. But if x is negative then we've mistakenly lost the negative sign.

With that in mind, here is a variant on your original query, with a form which retains the signum of a*omega*L outside the radical. Notice the broader equivalence, compared to your original target.

restart;

sol := (-v + sqrt(-4*a^2*R^2 + v^2))/(2*a*omega*L);

(1/2)*(-v+(-4*R^2*a^2+v^2)^(1/2))/(a*omega*L)

alt2 := -v/(2*a*omega*L) + sqrt((-4*R^2*a^2 + v^2)/(a^2*omega^2*L^2))/(2*signum(a*omega*L));

-(1/2)*v/(a*omega*L)+(1/2)*((-4*R^2*a^2+v^2)/(a^2*omega^2*L^2))^(1/2)/signum(a*omega*L)

combine(simplify(sol - alt2)) assuming real;

0

alt2 assuming a*omega*L>0;

-(1/2)*v/(a*omega*L)+(1/2)*((-4*R^2*a^2+v^2)/(a^2*omega^2*L^2))^(1/2)

alt2 assuming a*omega*L<0;

-(1/2)*v/(a*omega*L)-(1/2)*((-4*R^2*a^2+v^2)/(a^2*omega^2*L^2))^(1/2)

Download rad_comb_ex3.mw


Note: If, say, x=-3 then by convention sqrt(x^2) taken to be 3, not -3. That positive answer is known as the principal square root. See also the third paragraph here. Or see notes on the square-root case on principal value or n-th roots.

Is this the kind of thing you're trying to accomplish?

(I made no effort to optimize, or figure out what you're doing with tau at the end.)

I'll find a value for parameter beta that attains the goal g(10)=1. I'll do this in two ways, first using my own solving procedure and second using Doug Meade's Shoot package (see URL below).

restart

Digits:=15:

with(plots):

alias( U=u(x,y), V=v(x,y) ):

PDE:={ U*diff(U,x)+V*diff(U,y)-nu*diff(U,y$2)=0,

       diff(U,x)+diff(V,y)=0 }:

simsubs := eta(x, y) = y*sqrt((1/2)*u[0]/(nu*x)):

stream := psi(x, y) = sqrt(2*nu*x*u[0])*f(eta(x, y)):

Usubs := U = diff(rhs(stream), y):

Vsubs := V = -(diff(rhs(stream), x)):

ODE := simplify(subs(Usubs, Vsubs, simsubs, PDE)):

simsubs2 := solve(subs(eta(x, y) = eta, simsubs), {y}):

ODE := simplify(subs(simsubs2, ODE), symbolic):

FNS := {f(eta), g(eta), h(eta)};

{f(eta), g(eta), h(eta)}

ODE:={ diff(f(eta),eta)=g(eta),

       diff(g(eta),eta)=h(eta),

       diff(h(eta),eta)=-f(eta)*h(eta) }:

IC:={ f(0)=0, g(0)=0, h(0)=beta }:

BC:=g(10)=1:

 

findbeta := proc(b) local tt;
  if not b::numeric then return 'procname'(b); end if;
  tt := dsolve(eval(ODE union IC, beta=b),{f(eta),g(eta),h(eta)},numeric):
  eval(g(eta),tt(10));
end proc:

betaval := beta = fsolve(findbeta(beta)-1, beta=0..1);

beta = .469600090249286

Sf := dsolve(eval(ODE union IC, betaval),{f(eta),g(eta),h(eta)},numeric,output=listprocedure):


Check that g(10)=1 for this solution, computed using betaval.

eval(g(eta), Sf)(10);

HFloat(1.0000000000000007)

odeplot( Sf, [ [eta,f(eta)],

            [eta,g(eta)], [eta,h(eta)] ],0..4,
            labels=[`eta`,``], legend=["f","f '","f ''"],
            title="Blasius solution for flat-plate boundary layer" );

fp0:=subs( Sf, g(eta) ):

fp := proc(x)

  if not type(evalf(x),`numeric`) then

    'procname'(x);

  else

    fp0(x);

  end if;

end proc:

eta[`99%`] = fsolve(fp(eta) = .99, eta = 3 .. 4);

eta[`99%`] = 3.47188559058681

delta[`99%`] = subs(eta = rhs(%), combine(rhs(op(simsubs2))));

delta[`99%`] = 3.47188559058681*2^(1/2)/(u[0]/(nu*x))^(1/2)

tau := mu*simplify(diff(subs(simsubs, rhs(Usubs)), y));

(1/2)*mu*(nu*x*u[0])^(1/2)*((D@@2)(f))((1/2)*y*2^(1/2)*(u[0]/(nu*x))^(1/2))*2^(1/2)*u[0]/(nu*x)

tau[w] := subs(((D@@2)(f))(0) = (subs(Sf, h(eta)))(0), subs(y = 0, tau));

HFloat(0.234800045124643)*mu*(nu*x*u[0])^(1/2)*2^(1/2)*u[0]/(nu*x)


Now do the same thing using,
    http://www.math.sc.edu/~meade/maple/Shoot9/Shoot9.zip

libname:=cat(kernelopts(homedir),"/mapleprimes/Shoot/Shoot9"),libname:
with(Shoot);

[shoot]

infolevel[shoot] := 1:

S:=shoot( ODE, IC, BC, FNS, beta=1,

          abserr=Float(5,-7), output=listprocedure, method=taylorseries ):

shoot: Step #  1

shoot: Parameter values :  beta = 1
shoot: Step #  2
shoot: Parameter values :  beta = HFloat(0.4062401739572171)
shoot: Step #  3
shoot: Parameter values :  beta = HFloat(0.46805773973620646)

shoot: Step #  4
shoot: Parameter values :  beta = HFloat(0.4695991426488529)
shoot: Step #  5
shoot: Parameter values :  beta = HFloat(0.46959998836128514)


Check that g(10)=1 for this solution (approximately), computed using betaval.
(You could vary the options to shoot.)

eval(g(eta), S)(10);

.999999999972365

 

 

odeplot( S, [ [eta,f(eta)],

            [eta,g(eta)], [eta,h(eta)] ],0..4,
            labels=[`eta`,``], legend=["f","f '","f ''"],
            title="Blasius solution for flat-plate boundary layer" );

fp0:=subs( S, g(eta) ):

fp := proc(x)

  if not type(evalf(x),`numeric`) then

    'procname'(x);

  else

    fp0(x);

  end if;

end proc:

eta[`99%`] = fsolve(fp(eta) = .99, eta = 3 .. 4);

eta[`99%`] = 3.47188688127711

delta[`99%`] = subs(eta = rhs(%), combine(rhs(op(simsubs2))));

delta[`99%`] = 3.47188688127711*2^(1/2)/(u[0]/(nu*x))^(1/2)

tau := mu*simplify(diff(subs(simsubs, rhs(Usubs)), y));

(1/2)*mu*(nu*x*u[0])^(1/2)*((D@@2)(f))((1/2)*y*2^(1/2)*(u[0]/(nu*x))^(1/2))*2^(1/2)*u[0]/(nu*x)

tau[w] := subs(((D@@2)(f))(0) = (subs(S, h(eta)))(0), subs(y = 0, tau));

HFloat(0.23479999418064257)*mu*(nu*x*u[0])^(1/2)*2^(1/2)*u[0]/(nu*x)

NULL

Shoot_Blasius_solution_ac.mw

Is this the kind of effect you're after?

[edit:The solution(s) in terms of explicit radicals can be obtained by using the explicit option of the solve command. The simplification to your target form can be achieved via the evala command -- though there are a few longer incantations that also get it.]

restart

NULLeq1 := (R2+Rhi)*R1/(R1+R2+Rhi) = RloNULL

(R2+Rhi)*R1/(R1+R2+Rhi) = Rlo

 

NULLeq2 := R1*Rlo/(R1+Rlo)+R2 = RhiNULL

R1*Rlo/(R1+Rlo)+R2 = Rhi

evala({solve({eq1, eq2}, {R1, R2}, explicit)})

{{R1 = ((Rhi-Rlo)*Rhi)^(1/2)*Rlo/(Rhi-Rlo), R2 = ((Rhi-Rlo)*Rhi)^(1/2)}, {R1 = -((Rhi-Rlo)*Rhi)^(1/2)*Rlo/(Rhi-Rlo), R2 = -((Rhi-Rlo)*Rhi)^(1/2)}}


Download for_help_jrive_ac.mw

Maybe there is something nicer which "understands" the inert `%.` (as well as recognizes the distinct Matrix instances with equal entries as being equivalent), but it's less ad hoc.

ras5_v1_acc.mw

[edit] Here is very a minor revision to that, splitting to show the simplify step separately,

restart; Digits := 14: with(PDEtools): with(LinearAlgebra): with(plots):

pde:= diff(u(t,x),t$2) = diff(u(t,x),x$2):

 Subs:= {r1(t,x) = diff(u(t,x),t), r2(t,x)= diff(u(t,x),x)}:

eq1 := map(z->diff(z,t),<r1(t,x),r2(t,x)>) = Matrix(2,[0,1,1,0]) %.map(z->diff(z,x),<r1(t,x),r2(t,x)>):

#value(subs(Subs,eq1)):

GenerateStencil := proc(F,N,{orientation:=center,stepsize:=s,showorder:=true, showerror:=false}) loc...

eq2:=GenerateStencil(diff(r1(t,x),t),2,orientation=right,stepsize=s,showorder=false):
eq3:=GenerateStencil(diff(r2(t,x),t),2,orientation=right,stepsize=s,showorder=false):
eq4:=GenerateStencil(diff(r1(t,x),x),2,orientation=center,stepsize=h,showorder=false):
eq5:=GenerateStencil(diff(r2(t,x),x),2,orientation=center,stepsize=h,showorder=false):
 

eq6 := <rhs(eq2),rhs(eq3)> = Matrix(2,[0,1,1,0]) %. <rhs(eq4),rhs(eq5)>:

eq7 := [seq(seq(r1(t+i*s,x+j*h)= r1[n+i,m+j],i=0..1),j=-1..1)]:
eq8 :=[seq(seq(r2(t+i*s,x+j*h)= r2 [n+i,m+j],i=0..1),j=-1..1)]:

eq9 := map(z->subs(eq7,eq8,z),eq6):

eq10 := [seq(seq(r1(t+i*s,x+j*h)= r1[n+i,m+j]-E1[n+i,m+j],i=0..1),j=-1..1)]:

eq11 :=[seq(seq(r2(t+i*s,x+j*h)= r2 [n+i,m+j]-E2[n+i,m+j],i=0..1),j=-1..1)]:

eq12 := map(z->subs(eq10,eq11,z),eq6):

eq13 := simplify(eq12-eq9,'mindeg');

(Vector(2, {(1) = (E1[n, m]-E1[n+1, m])/s, (2) = (E2[n, m]-E2[n+1, m])/s})) = `%.`(Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 1, (2, 2) = 0}), Vector(2, {(1) = (1/2)*(-r1[n, m-1]+E1[n, m-1]+r1[n, m+1]-E1[n, m+1])/h, (2) = (1/2)*(-r2[n, m-1]+E2[n, m-1]+r2[n, m+1]-E2[n, m+1])/h}))-`%.`(Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 1, (2, 2) = 0}), Vector(2, {(1) = (1/2)*(-r1[n, m-1]+r1[n, m+1])/h, (2) = (1/2)*(-r2[n, m-1]+r2[n, m+1])/h}))

factor(subs[eval](`%.`=`*`,subsindets(subsindets(eq13,Matrix,
                                                 u->freeze(convert(u,listlist))),
                                      Vector,u->freeze(convert(u,list))))):
eq14 := subsindets(subsindets(subsindets(thaw(%),And(`*`,satisfies(u->ormap(type,u,list))),
                                         `%.`@op),listlist,Matrix),list,Vector);

(Vector(2, {(1) = (E1[n, m]-E1[n+1, m])/s, (2) = (E2[n, m]-E2[n+1, m])/s})) = `%.`(Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 1, (2, 2) = 0}), Vector(2, {(1) = -(1/2)*(-r1[n, m-1]+r1[n, m+1])/h+(1/2)*(-r1[n, m-1]+E1[n, m-1]+r1[n, m+1]-E1[n, m+1])/h, (2) = -(1/2)*(-r2[n, m-1]+r2[n, m+1])/h+(1/2)*(-r2[n, m-1]+E2[n, m-1]+r2[n, m+1]-E2[n, m+1])/h}))

simplify(%);

(Vector(2, {(1) = (E1[n, m]-E1[n+1, m])/s, (2) = (E2[n, m]-E2[n+1, m])/s})) = `%.`(Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 1, (2, 2) = 0}), Vector(2, {(1) = (1/2)*(E1[n, m-1]-E1[n, m+1])/h, (2) = (1/2)*(E2[n, m-1]-E2[n, m+1])/h}))

Download ras5_v1_accc.mw

The technique might be adapted to handle several kinds of more involved Matrix/vector expressions with inert `%.`. That's related to a major reason why I made the code so involved -- for generality.

If results is a list, set, Vector/Array, or expression sequence (of, say, 7 elements) then you could do any of,

    return  results[2];

    return results[3..6];

    return results[2],results[4],results[5];

    return results[1],results[7];

Here are a few edits. I've added comments as to why I've made the changes.

Note that I've also allowed your utility formatters (Code-Edit region) to continue to function like the original.

(In the actual Maple GUI the units render more nicely. The double-braces appearing around units below are an artefact of this website's old backend.)

restart

kernelopts(version)

`Maple 2024.0, X86 64 LINUX, Mar 01 2024, Build ID 1794891`

conv2ftin:=proc(arg)

NULL

 

NULLNULL

totvol := proc (ht) options operator, arrow; Pi*.4^2*ht^2*ht+(1/3)*Pi*.4^2*ht^2*((2/3)*ht) end proc

proc (ht) options operator, arrow; .6143558967*ht^3 end proc

The procedures defined before this package is loaded
(ie. your utilities, as well as `totvol`) get the global
bindings for arithmetic.
That allows the utilities to work as intended, and
it allows the body of `totvol` to automatically simplify
at creation time, as well as to print as you wanted.

with(Units)

Automatically loading the Units[Simple] subpackage
 

Utilize the foot-pound-second system instead of SI.

UseSystem(FPS)

totvol(2.369133117)

8.169367285

fsolve(totvol(ht) = 8.169367284, ht)

2.369133117

volReq := 130*Unit('gallon'['US_liquid'])

130*Units:-Unit(gal)

volReq := 1.2*volReq

156.0*Units:-Unit(gal)

Alternate to next pair of commands below.
This produces result in ft, without having to first
convert volReq from gal to ft^3. For the FPS system
length dimensions computations simplify to ft.

volReq := convert(volReq, 'units', ft^3)

20.85416667*Units:-Unit(ft^3)

By default fsolve computes real roots.
It's suitable for your floating-point example.
Note that the result is in ft, not meters.

Also, by assigning directly to ht__main_tank we
avoid assigning to name `ht` which is used as
dummy name in a few places. This allows these
parts of your worksheet to be recomputed out-of-order,
since name ht is not clobbered.

ht__main_tank := fsolve(totvol(ht) = volReq, ht)

3.237856544*Units:-Unit(ft)

 

With the computed height in ft, we don't need
to convert all of these results separately.

ht__pan := (2/3)*ht__main_tank; rad := .4*ht__main_tank; total_height := ht__main_tank+ht__pan

2.158571029*Units:-Unit(ft)

1.295142618*Units:-Unit(ft)

5.396427573*Units:-Unit(ft)

Note that your utilities still function as intended.
(They don't continue to work properly if you just follow the
Answers suggesting that all you need to do it utilize `unapply`...)

conv2ftin(ht__main_tank); conv2ftin(ht__pan); conv2ftin(rad); conv2ftin(total_height)

3.*Units:-Unit(ft)+2.854278528*Units:-Unit(`in`)

2.*Units:-Unit(ft)+1.902852348*Units:-Unit(`in`)

1.*Units:-Unit(ft)+3.541711416*Units:-Unit(`in`)

5.*Units:-Unit(ft)+4.757130900*Units:-Unit(`in`)

convert(totvol(ht__main_tank), 'units', gallons)

156.0000001*Units:-Unit(gal)

Download Tank_Design_Calculation_-_Units_Question_(v00)_ac.mw

You can utilize the assuming facility to place separate assumptions during individual computations.

That allows you to compute under the assumption m::integer, but with no assumption on alpha.

For example, using Maple 2022 because the OP did,

restart;

kernelopts(version);

`Maple 2022.2, X86 64 LINUX, Oct 23 2022, Build ID 1657361`

expr2 := 4*cos(Pi*m+2*alpha*Pi)/(Pi*(2*m + 1));

4*cos(2*Pi*alpha+Pi*m)/(Pi*(2*m+1))

combine(expand(expr2)) assuming m::integer;

4*(-1)^m*cos(2*alpha*Pi)/(Pi*(2*m+1))

Download assuming_trig_ex.mw

It sounds as if you want nice pretty-printed ways to both enter and print those (2D Input and Output), and not just the means to compute.

For a more involved attempt you could look at this old Answer thread (and also some later revisions in worksheets attached to replies in the thread.)

For a preview, here's some simple functionality along those lines, again with 2D Input/Output.

restart;
`convert/degpolar`:=proc(x)
   degpolar(abs(x),180/Pi*argument(x));
end proc:
`print/degpolar`:=proc(ma,ar)
   uses Typesetting;
   mrow(Typeset(`&angle;`(ma,ar)),mi("&deg;"));
end proc:
`evalc/degpolar`:=proc(ma,ar)
   ma*(cos(ar*Pi/180)+I*sin(ar*Pi/180));
end proc:

`&angle;` := degpolar

degpolar

 

 

 

a := 1/2+sqrt(3)/2*I;

1/2+((1/2)*I)*3^(1/2)

ans := convert(a, degpolar);

degpolar(1, 60)

evalc(ans);

1/2+((1/2)*I)*3^(1/2)

 

For the next I've used the symbol from the Operators Palette.

 

q := `&angle;`(230, 0)/`&angle;`(92, -53.13)

degpolar(230, 0)/degpolar(92, -53.13)

convert(evalc(q), degpolar)

degpolar(2.500000001, 53.12999999)

Download degpolar_ch.mw

You specifically inquired about implicitplot3d, so in coding below I will address that in particular. It's documented [123] and not too difficult to color/shade a surface constructed by plot3d, etc. Hence if you could break up your implicit x-y-z region into explicitly defined portions then you could use another 3d plotting command. So I'm supposing that you wouldn't be asking about implicitplot3d if you had only an example that was subject to such recourse.

I'm just going to address coloring (arbitrarily) with implicitplot3d, since dharr's already treated your particular surface examples with a split into explicitly formulated portions.

You can in fact utilize a float[8] Array to color the triangles by which the GUI renders the surface, when it displays an ISOSURFACE plotting structure generated by the implicitplot3d command. You can use such an Array to specify coloring by either RGB or HSV values.

But the big problem is that those triangles don't match the set of x-y-z points stored as data in the ISOSURFACE. The ISOSURFACE values are just for the regularly spaced points in the ranges. The location of the vertices of those triangles are a result of the interpolating process by which the GUI's plot driver figures out where the surface might be, and we are not privy to that. So we lack programmatic access to the x-y-z values for the vertices of the GUI-rendered triangles.

Even the number of the mesh points of stored function values doesn't match the number of triangles in the rendered surface. And that number of triangles is also a function of the proportion of the plotting window which the surface occupies.

Here's an example. This code does not form those triangles. The GUI does that. This code merely adds a COLOR substructure to a usual implicitplot3d result, to specifiy hue values for its rendered triangles.

subsindets(plots:-implicitplot3d(x^2+y^2+z^2=3.9,x=-2..2,y=-2..2,z=-2..2,grid=[6,6,6],
                                 lightmodel=none,orientation=[60,75,0]),
           specfunc(ISOSURFACE),
           Iso->ISOSURFACE(op(Iso),
                           COLOUR(HSV,Array(1..1528,1..3,
                                            (ijk,hsv)->
                                              `if`(hsv=1,
                                                   piecewise[1.0](irem(ijk,3)=0,0.1,
                                                                  irem(ijk,3)=1,0.4,
                                                                  irem(ijk,3)=2,0.7),
                                                   1.0),
                                            datatype=float[8]))));

 

 

 

Download ISO_col_0.mw

note: If you changed that to utilize grid=[6,6,7] instead then in Maple 2016.2 and later you could see an example of the GUI message spit out when the coloring Array doesn't have the right number of entries.

As it happens that COLOR substructure's Array can match either the number of segments of the number of vertices (whichever is closest, I suspect). And the GUI spits out a message when there's a bad mismatch. But in neither case do we know the x-y-z coordinates of those vertices, and so we don't know how to populate that COLOR structure's Array. (If I had such a means I'd ditch ISOSURFACE myself, and do the math and draw those portions by my own construction, with full control...)

If I had to distinguish such from a mere procedure call then I'd refer to such as an indexed procedure call.

Now, an index on the base name of a procedure call only would only have a special effect if someone had actually written the original Library procedure to make use of such an index. Otherwise it'd act as it would without the indexing, ie. the index would be simply ignored.

Run your combine and simplify examples involving indexed procedure calls again, but now without the index, and you should get the same results as with the index. Those routines are not coded to recognize and act according to such an index.

Eg, the index [trig] on combine is not what makes this example work as it does. It's just the same result as if no optional arguments were supplied.

combine[trig](exp(sin(a)*cos(b))*exp(cos(a)*sin(b)));

           exp(sin(a + b))

combine(exp(sin(a)*cos(b))*exp(cos(a)*sin(b)));

           exp(sin(a + b))

I don't really understand why you thought that all your modified examples of indexed procedure calls would do something special. The documentation Calling Sequence, Description, and Examples don't indicate that they would.

As for your examples with extra passed options, perhaps additional commentary or examples might make it more clear on the Help page ?combine that the choice of such supplied options may restrict (by virtue of specifying) what manner of combining gets done. Notice how in the following examples the results are each combined only in the sense of the specified option,

combine(exp(sin(a)*cos(b))*exp(cos(a)*sin(b)),[exp]);

combine(exp(sin(a)*cos(b))*exp(cos(a)*sin(b)),[trig]);

combine(exp(sin(a)*cos(b))*exp(cos(a)*sin(b))
        +exp(sin(a)*cos(b)+cos(a)*sin(b)),[trig]);

As for your evalf[n] mention, evalf has been specially written to take account of an index on the called procedure name, and use it like an option. That is documented.

You could utilize the known option of dsolve, using a wrapper around Mt that is just a procedure (that returned unevaluated for nonnumeric input).

The known option is a usual way to inform dsolve that you have a helper (black box) function.

The dsolve command (with its known option) does not appear to recognize directly the kind of appliable object that Interpolation returns. Hence my kludge with a constructed wrapping operator/procedure.

You can compare the result with simply using a linear expression as the rhs driver.

restart;

M := Array(1 .. 10, 1 .. 2):
for i to 10 do
    M[i, 1] := i;
    M[i, 2] := 3*i;
end do:

Mt := Interpolation:-LinearInterpolation(M):

E := unapply('Mt'(t),t,numeric):

diffeq := D(C)(t) = E(t);

(D(C))(t) = E(t)

solA := dsolve({diffeq, C(0) = 0}, {C(t)}, numeric, known=E):

plots:-odeplot(solA,[t,C(t)]);

diffeqB := D(C)(t) = 3*t;
solB := dsolve({diffeqB, C(0) = 0}, {C(t)}, numeric):

(D(C))(t) = 3*t

plots:-odeplot(solB,[t,C(t)]);


Download ds_kn.mw

A variety of ways, with fun still defined outside the procedure.

restart;

fun := x^2+y^2;

x^2+y^2

xt := -1: yt := 2:


Your procedure could evaluate the expression fun by
replacing the global x,y (in fun) by the values of the
procedure's parameters x,y.

Also, more efficiently doing just a single call to the procedure
and making a multiple assignment of its returned sequence.

temp_proc := proc(x, y)
  local lfun, out1, out2, out3;
  lfun := eval(fun,[:-x=x,:-y=y]);
  if x > 0 then
    out1,out2,out3 := lfun, 2*lfun, k*lfun;
  elif x <= 0 then
    out1,out2,out3 := lfun, -2*lfun, -k*lfun;
  end if;
  return out1, out2, out3;
end proc:

out1_fin,out2_fin,out3_fin := temp_proc(xt, yt);

5, -10, -5*k


The same as before, but assigning the sequence of three results
to a single local out.

temp_procB := proc(x, y)
  local lfun, out;
  lfun := eval(fun,[:-x=x,:-y=y]);
  if x > 0 then
    out := lfun, 2*lfun, k*lfun;
  elif x <= 0 then
    out := lfun, -2*lfun, -k*lfun;
  end if;
  return out;
end proc:

out1_fin,out2_fin,out3_fin := temp_procB(xt, yt);

5, -10, -5*k


The same as the last, but without assigning to a
local out. The last result in the proc is what gets returned.

temp_procC := proc(x, y)
  local lfun;
  lfun := eval(fun,[:-x=x,:-y=y]);
  if x > 0 then
    lfun, 2*lfun, k*lfun;
  elif x <= 0 then
    lfun, -2*lfun, -k*lfun;
  end if;
end proc:

out1_fin,out2_fin,out3_fin := temp_procC(xt, yt);

5, -10, -5*k


The same a the last, but using a passed operator (formed
just once from the expression fun) instead of evaluating
to replace the global name instances.

Note that Fun is only called once, inside the procedure,
which is efficient.

 

One point that might be key to you is that you might already
have fun as an expression computed elsewhere. The unapply constructs
the operator from that.

 

temp_procD := proc(x, y, F)
  local lfun;
  lfun := F(x,y);
  if x > 0 then
    lfun, 2*lfun, k*lfun;
  elif x <= 0 then
    lfun, -2*lfun, -k*lfun;
  end if;
end proc:

Fun := unapply(fun,x,y);

proc (x, y) options operator, arrow; x^2+y^2 end proc

out1_fin,out2_fin,out3_fin := temp_procD(xt, yt, Fun);

5, -10, -5*k

 

Like the previous, but without passing in the operator.

 

temp_procE := proc(x, y)
  local lfun;
  lfun := Fun(x,y);
  if x > 0 then
    lfun, 2*lfun, k*lfun;
  elif x <= 0 then
    lfun, -2*lfun, -k*lfun;
  end if;
end proc:

Fun := unapply(fun,x,y);

proc (x, y) options operator, arrow; x^2+y^2 end proc

out1_fin,out2_fin,out3_fin := temp_procE(xt, yt);

5, -10, -5*k

 

Download fun_var.mw

f := x -> (1-k*x)/(1+x^2):

Student:-Calculus1:-Tangent(f(x), x=3);

(2/25)*k*x-(3/50)*x-(27/50)*k+7/25

 

And, optionally,

 

collect(%, x);

x*((2/25)*k-3/50)-(27/50)*k+7/25

Download Tg_ex.mw

You didn't state how wide you wanted the columns.

The following shows that your stated goal of getting the columns centered at 1,2,3 can be done using only one extra (option) value chosen by you. Here the other options can be just functions of that -- see last example.

You could use 0 < width <= 1 below.

restart;

x:= Vector([1,2,3]):

Statistics:-ColumnGraph(x,'offset'=0.5,'distance'=0,'width'=1);

Statistics:-ColumnGraph(x,'offset'=1-0.7/2,'distance'=1-0.7,'width'=0.7);

Download CG_ex.mw

2 3 4 5 6 7 8 Last Page 4 of 313