acer

32303 Reputation

29 Badges

19 years, 307 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@janhardo I put it back in as the original float, to avoid confusion.

With our hands on the denominator (unexpanded, and also the "common" term), then we could handle this particular example like so.

Zin := ((Rsp + (omega*L2 - omega02^2*L2/omega)*I)
       *(Rpp + (omega*L1 - omega01^2*L1/omega)*I) + 0.01*omega^2*L1*L2)
       /(Rsp + (omega*L2 - omega02^2*L2/omega)*I);

((Rsp+I*(omega*L2-omega02^2*L2/omega))*(Rpp+I*(omega*L1-omega01^2*L1/omega))+0.1e-1*omega^2*L1*L2)/(Rsp+I*(omega*L2-omega02^2*L2/omega))

That denominator, untouched and not expanded...

FF := 1/select(type,Zin,anything^(negint));

Rsp+I*(omega*L2-omega02^2*L2/omega)

map(`/`,Zin*FF,FF);

Rpp+I*(omega*L1-omega01^2*L1/omega)+0.1e-1*omega^2*L1*L2/(Rsp+I*(omega*L2-omega02^2*L2/omega))

Download collect_ex6.mw

I ought to point out that this does not systematically check for the particular form that allows this to work. Without such a systematic process I label this as ad hoc.

Since the cited target result is not the simplest form possible then defining such a systematic process ought to be preceded by a clear definition of the candidate form. (I think I see one.)

On the other hand, the choice of variables that were used earlier is also an ad hoc aspect.

@nm 

Yes, Mma's simplification is ok. But the following is even simpler.

Rpp+(I*(omega-omega01^2/omega)+0.1e-1*omega^2*L2/(Rsp+I*(omega-omega02^2/omega)*L2))*L1


I can get that with a minor variant. (I didn't before; it was not the OP's target.)


ps. quibble: I notice that you changed the 0.01 to 100, in your Mma example.

@jrive I'm sorry that my effort is so looked-down upon.

Personally, I find the following to be slightly complicated but not advanced.

FF := 1/select(type,Zin,anything^(-1)):
simplify(Zin,{FF=freeze(FF)}):
thaw(collect(%,[Rpp,L2],simplify));

1.000000000*Rpp+0.1e-1*omega^2*L1*L2/(Rsp+I*(omega*L2-omega02^2*L2/omega))+I*L1*(omega^2-1.*omega01^2)/omega

or,

FF:=1/select(type,Zin,anything^(-1)):
simplify(Zin,{FF=FFF}):
subs(II=I,FFF=FF,collect(subs(I=II,%),[Rpp,II,omega]));

1.000000000*Rpp+I*(1.000000000*omega*L1-1.000000000*omega01^2*L1/omega)+0.1000000000e-1*omega^2*L1*L2/(Rsp+I*(omega*L2-omega02^2*L2/omega))

I made an effort to not have those irritating float artefacts.

I also made an effort to explain the extra steps I took to work around issues (eg. collecting wrt I, the float contagion, etc).

On the bright side at least it isn't one of those very poor and ad hoc solutions involving op and by-hand/eye picking of operand numbers/positions. That's a nearly useless approach!

You tagged this with simplify. But if you're attempting some notion of simplification then why does your target form not have the L1 (or L2) factored out of some subexpressions, eg,
   omega*L1 - omega01^2*L1/omega
instead of,
    L1*(omega - omega01^2/omega)
?

I'm just curious. The target can be attained using collect as the main tool. But I want to be sure that I am using the right target.

Also, it looks as if the Rac term was just some typo.

@C_R The particular error message shown by the OP is not the primary issue. The outer integration won't succeed as intended without proper syntax, which it lacked by having an fdiff call as first argument and a bare range as second argument. In M2024.1 one might get back, say, an unevaluated int call, say, but such a result is not meaningful or useful.

Numeric integration is wanted in this case, as part of this described exercise. And one can enable that here by one of these 1) using the numeric option to lowercase int, 2) making the upper bounds passed to int be numeric floats, or 3) invoking evalf(Int(...)) instead. So that's one aspect.

Another aspect is to also repair the syntax for the arguments passed, and fixing that is required if a float result is wanted. That can be achieved for all three of the above approaches, by using an operator for the first argument so that the x value used in each fdiff call is an actual numeric value, as I mentioned in my Answer. That happens to succeed (slowly), even with stock options, for this given example and time-value t=0.

Ie, using,
   x->fdiff(sol(x, t), [t = 0]), a .. b
as first arguments for all three,

kernelopts(version);

`Maple 2024.1, X86 64 LINUX, Jun 25 2024, Build ID 1835466`

Using option numeric to the int call.
restart;
with(PDEtools): with(plots): with(DEtools):
pde := diff(v(x, t), t, t) = diff(v(x, t), x, x):
f := xi -> exp(-xi^2):
a := -10: b := 10: dx := 1/50: t_final := 10:
pds := pdsolve(pde, {v(a, t) = 0, v(b, t) = 0,
               v(x, 0) = f(x + 5), D[2](v)(x, 0) = -eval(diff(f(x), x), x = x + 5)},
               numeric, range = a .. b, time = t, spacestep = dx):
sol_proc := rhs(pds:-value(output = listprocedure)[3]):
sol := (x, t) -> piecewise(a < x and x < b, sol_proc(x, t), 0):
int(x->fdiff(sol(x, t), [t = 0]), a .. b, numeric);

0.2324739513e-10

Using float values in the range in the int call.
restart;
with(PDEtools): with(plots): with(DEtools):
pde := diff(v(x, t), t, t) = diff(v(x, t), x, x):
f := xi -> exp(-xi^2):
a := -10.0: b := 10.0: dx := 1/50: t_final := 10:
pds := pdsolve(pde, {v(a, t) = 0, v(b, t) = 0,
               v(x, 0) = f(x + 5), D[2](v)(x, 0) = -eval(diff(f(x), x), x = x + 5)},
               numeric, range = a .. b, time = t, spacestep = dx):
sol_proc := rhs(pds:-value(output = listprocedure)[3]):
sol := (x, t) -> piecewise(a < x and x < b, sol_proc(x, t), 0):
int(x->fdiff(sol(x, t), [t = 0]), a .. b);

0.2324739513e-10

Using evalf(Int(...)) instead of int(...).
restart;
with(PDEtools): with(plots): with(DEtools):
pde := diff(v(x, t), t, t) = diff(v(x, t), x, x):
f := xi -> exp(-xi^2):
a := -10: b := 10: dx := 1/50: t_final := 10:
pds := pdsolve(pde, {v(a, t) = 0, v(b, t) = 0,
               v(x, 0) = f(x + 5), D[2](v)(x, 0) = -eval(diff(f(x), x), x = x + 5)},
               numeric, range = a .. b, time = t, spacestep = dx):
sol_proc := rhs(pds:-value(output = listprocedure)[3]):
sol := (x, t) -> piecewise(a < x and x < b, sol_proc(x, t), 0):
evalf(Int(x->fdiff(sol(x, t), [t = 0]), a .. b));

0.2324739513e-10


Download pde_num_e2.mw

@mmcdara Thanks, that's the kind of thing I'd surmised. But the numeric requirements are imprecise.

My concern is that such a process is not a test of only the numeric pde solver. Its results rely upon not only behavior of the three different & separate numeric solvers I mentioned but also upon their interplay.

Failure to attain some (as yet unstated) degree of numeric accuracy cannot therefore be a proper gauge of numeric performance of any of sole component part of that process. Eg. failure might say nothing about the numeric pde solver alone.

I supposed that accuracy of the integration step would need be coarser than that of the diffentiation step, that in turn being coarser than that of the pde scheme. But the particulars of that might need be problem specific. And attaining an upfront, explicit, accuracy goal might need per-problem tuning at any/all of those levels.

@GFY By the way, it is also possible to do it with fsolve (as a modification/correction of the earlier attempt). Here, I augment with undefined, for the absent root values.

You can also scale any you choose, eg. by 3.

I use option remember to avoid calling fsolve multiple times for the same input (ie. same omega__v, or r, value).

You could make this more sophisticated.

restart;

plots:-setoptions(size=[500,200]);

det_num := -0.9633005975*Omega^6 + (0.9630831768*omega__v^2 + 0.9337694145)*Omega^4 + (-0.9338923707*omega__v^2 - 0.01494064902)*Omega^2 + 0.01494064902*omega__v^2 = 0;

-.9633005975*Omega^6+(.9630831768*omega__v^2+.9337694145)*Omega^4+(-.9338923707*omega__v^2-0.1494064902e-1)*Omega^2+0.1494064902e-1*omega__v^2 = 0

ll := proc(varr)
  local i, res; option remember;
  res := [fsolve(eval(det_num, omega__v=varr), Omega=0..infinity)];
  return [res[], seq(undefined, i=1..3-nops(res))];
end proc:

ll(0.8);

[.1275587179, .8003574462, .9758873053]

plots:-display(
  plot(x->ll(x)[1],0..2),
  plot(x->ll(x)[2],0..2),
  plot(x->ll(x)[3],0..2));

det_num := 0.7801238870*Omega^4 + (-9.794000000/r - 1340.585417)*Omega^2 + 13129.69357/r = 0

.7801238870*Omega^4+(-9.794000000/r-1340.585417)*Omega^2+13129.69357/r = 0

ll2 := proc(varr)
  local i, res; option remember;
  res := [fsolve(eval(det_num, r=varr), Omega=0..infinity)];
  return [res[], seq(undefined, i=1..2-nops(res))];
end proc:

plots:-display(
  plot(x->ll2(x)[2],0.01 .. 0.1),
  plot(x->3*ll2(x)[1],0.01 .. 0.1));

 

Download fuxian_ac2.mw

@GFY Firstly, note that there is a range/set of values of omega_v for which there is only one positive root Omega of det_num. That is why there is a middle section with no upper root. It's not that the method fails to properly compute them; they don't exist in that region.

[seq(RootOf(det_num, Omega, index = i), i = 2 .. 3)];
plot(%, omega__v = 0.95 .. 1.01, size = [400, 400]);

Secondly, you are not being clear about what you mean by the term "root". Even if you were able to construct explicit formulas (with radicals) those need not match up with the apparent "curves" in the way you expect. (That can happen even for a parametrized quadratic.) The numeric rootfinding method in play here involves generation of all the roots (per omega__v value), which you may sort. That's what you had originally; you computed all roots (for each given omega__v) and then you applied sort. And that's essentially what evalf/RootOf does here too.

It's a mistake to think that I did something fundamentally different than what you were trying earlier. Using indexed RootOf was just easier than repairing your procedure to inject NULL in the cases/region of less actual positive roots.

You might be able to do numeric rootfinding by a homotopy (or other) approach instead, which gets curves of solid color. If you want color to be preserved across some crossing point of roots then you should at least define that concept.

Another way, with uniform color,

plots:-implicitplot(det_num,omega__v=0.0..2.0,Omega=0..2,
                    rangeasview, gridrefine=4, signchange=false);

You might be able track a "root" curve through a simple crossing/switch by matching derivatives.

@maple2015 You can alter that line to,

   newintg :=  collect(intg,lhs~(inter),uu->simplify(uu,size))

so that lhs~(inter) becomes the list [r,t] or [u,v], etc, as appropriate.

ps. In your original worksheet case 3 had a different substitution, and the numeric integration did not complete successfully on my machine. If you have a revision, you could attach the actual worksheet here.

I have deleted a duplicate Question thread on this.

If you have followup queries on this, then please put them here instead of spawning wholly separate new Question threads.

@Carl Love It's available by doing right-click on his RUN button and choosing the item "Edit Click Action".

Part of the time is spent in a numeric integration in variables r & t. It seems faster to do,
  newintg := simplify(intg,size):
nd then,
  evalf(Int(newintg,inter,digits=15,epsilon=1e-7,method=_cuhre))
without setting Digits=7 (it's Digits=30 earlier, though I don't know exactly why).

It may also be quicker to do,
   newintg := collect(intg,[r,t],uu->simplify(uu,size)):
than to hit the whole with simplify(...,size).

The other part of the cost (about half, maybe) seems to be initial construction of formulas, in particular the line,
   st := [seq(coeftayl(Eq, y = 0, j), j = 0 .. M-2)]
but I haven't looked at that.

@Susana30 You still haven't explained what 2D plot you want.

Maybe it's one of these...
some_plots.mw

(check for mistakes. I don't really know what ranges you want, or if any restricted.)

There are several different ways to construct these (or project them). I duplicated some computation, to illustrate.

@Andiguys 

Your worksheet has s and s_1 as tables. And you know the tau0 values used. (You can also access the indices of the tables programmatically.) So constructing either a mx2 Matrix or a list-of-lists should be straightforward, and you can pass such into the plot command.

Eg,

plot([seq([uu, s[uu][1]], uu = sort([indices(s, nolist)]))])

In my Maple 2019.2 the optimal values for s_1 don't really vary (except for some cruft artefacts of the HFloats displayed). I doubt that it's working as you intend. Increasing Digits:=30 at the very start doesn't help much. It's not clear to me whether tau0 is actually playing the actual role in your objectives that you want.

I'd suggest that the optimization itself might need fixing here. But it lacks comments and proper explanation.


   for tau0 in {0.3, 0.35, 0.4, 0.45, 0.5, 0.6} do
     ...
   end do

First 27 28 29 30 31 32 33 Last Page 29 of 591