acer

32348 Reputation

29 Badges

19 years, 329 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@tomleslie Various approaches, including Lagrange form interpolating polynomial, series solution (about x=0) via IVP reformulation, and Chebychev-Pade rational polynomial approximation.

restart;

de := sin(1)*(diff(y(x), x$2))+(1+cos(1)*x^2)*y(x) = -1:
cond := y(-1) = 0, y(1) = 0:

Digits:=15;
plist:=dsolve({cond, de}, y(x), numeric, output=listprocedure):
me:=dsolve({cond, de}, y(x), numeric, output=mesh)[2,1];

Digits := 15

 

Matrix(%id = 18446884398932593470)

(1)

 

CFy := CurveFitting:-PolynomialInterpolation(me[..,1..2],
                                             x,form=Lagrange):
degree(CFy,x);

14

(2)

P:=Student:-NumericalAnalysis:-PolynomialInterpolation(convert(me[..,1..2],
                                                               listlist),
                                                       method=lagrange):

NAy:=Student:-NumericalAnalysis:-Interpolant(P,independentvar=x):
degree(NAy,x);

14

(3)

sort(fnormal(expand(CFy-NAy)));

0.

(4)

Y:=eval(y(x),plist):

plot([abs(Y(x)-CFy)], x=-1..1, size=[600,200], view=0..1e-8);

 

plot([abs(Y(x)-NAy)], x=-1..1, size=[600,200], view=0..1e-8);

 

TLsol:=dsolve({de,cond}, numeric):
Order:=15:
TLy:=sort(evalf(rhs(convert(
             dsolve({de,
                     y(0)=rhs(TLsol(0)[2]),
                     D(y)(0)=rhs(TLsol(0)[3])},
                     y(x), 'series'), polynom)))):
degree(TLy,x);

14

(5)

plot([abs(Y(x)-TLy)], x=-1..1, size=[600,200], view=0..1e-8);

 

CPy:=eval(numapprox[chebpade](eval(y(x),plist)(x),
                              x=-1..1,[14,1]),T=orthopoly[T]):
degree(CPy,x);

14

(6)

plot([abs(Y(x)-CPy)], x=-1..1, size=[600,200], view=0..1e-8);

 

 

 


Download bvp_interp.mw

@ligonberry 

I don't recall an earlier question where someone asked for the explicit interpolant for a BVP.

I should ask why you want to get such a thing, assuming I'm interpreting you rightly. Is it for export, to be used in another program? If not then why do you want it?

There would be few reasons to want it if you're just going to be operating in Maple itself. The procedure returned by dsolve(...,numeric) is designed to utilize the computed mesh of interpolation points efficiently, on demand when given specific numeric input for the independent variable.

How would you feel about, say, an explicit symbolic interpolant constructed from the data returned by supplying the output=mesh option for dsolve(...,numeric)?

IIRC the usual procedure returned by dsolve(...numeric) for BVPs will itself call internals like `dsolve/numeric/lagrange` or `dsolve/numeric/hermite`, to efficiently generate the solution at a queried numeric point. Perhaps we could construct the explicit interpolant based on this usage.

@ligonberry Mariusz gave you Maple code which solves it numerically, as the first 5 lines of dark text in his Answer.

@Adam Ledger You are wrong when you claim,

   "Also because it is a simplification that is concealed in an inner level of
    abstraction of evalf, i need to know the showstat for evalf so that I can
    further understand the causation of the invalid computation"

The exact result from simplify is not correct. It's already wrong before it gets into evalf. You don't need to see the inner workings of evalf to figure out the cause of the discrepency. It's a branch-cut problem related to simplify, not a problem in computing to floating-point under evalf.

Here, simplify is behaving as if it's not taking into account that LambertW(1,Z) might not be real.

@Adam Ledger I shortened the example to try and illustrate more succinctly that the difference is due to an invalid simplification.

@SGJanssens Is there an attachment?

I'm not against a Post, if you have results to share. (You could also Branch off a Post as well.)

Please stop converting this from Question to Post, whoever is doing so.

G := tanh(h)^(1/2)*(exp(4*h)-1)^(1/2)*exp(-h)+exp(h)-exp(-h):

solve( G, [h] );   # ouch
                      [[h = h]]

@mmcdara 

Here are a couple of things that I sometimes do when using subsindets, as ways to check in advance whether the type (2nd argument) will target the subexpressions that I'm interested in modifying.

1) Before I call   subsindets(expr, sometype, myaction)  I may first try   indets(expr, sometype)

2) Before I call   subsindets(expr, sometype, myaction)  I may first try   subsindets(expr, sometype, K) where K is just an unassigned name.

BTW, I think that it should be not very hard to programmatically replace the ellipses with rectangles (of just one size, but corrected to enclose the longest vertex label at a given font size), in the result from DrawGraph. Trickier, but likely still possible, would be making each rectangle be sized to match the corresponding vertex label. Unfortunately I don't have time to try that these days.

@gkokovidis Clever of you to figure out the goal.

And, sure, it can also be done with Vectors (as I'm sure you know).

restart:
x:=<seq(0..20,0.5)>:
y:=map(i->5*i^2+3, x):
xy := <x|y>:
plot(xy);

@ALIKHADEMI Please show us an explicit example of what you're trying to do. Your description is too vague.

@mmcdara Even in Maple 2015.2 one can pass the option linestyle=dash straight to the DrawGraph command.

Doing so returns a PLOT structure with LINESTYLE as a global substructure, ie. not specific to any particular POLYGONS substructure.

And it's still possible to use the handy subsindets command to manipulate the individual POLYGONS substructures.  (In this case the edges are POLYGONS with a list of two lists as the first operand, and in Maple 2015 the vertex boxes are POLYGONS with a list of more lists as first operand).  So an appropriate type as second argument to subsindets allows finer targeting.

Graphfun_2015.mw

@Kitonum It may sometimes be possible to condition the first argument to solve, using discont.

These two simple examples don't produce a RealRange from discont, but in such an event the points and ranges could be assembled slightly differently.

restart;
f:=ln(x);
                                 f := ln(x)

g:=evalc(Im(f)):
discont(f,x);
                                     {0}

P:=seq(Or(x>p,x<p), p=discont(f,x));
                            P := Or(0 < x, x < 0)

sols:=[solve( {g, P}, {x} )];
                              sols := [{0 < x}]

restart;
f:=ln(x)*ln(x-3);
                            f := ln(x) ln(x - 3)

g:=evalc(Im(f)):
discont(f,x);
                                   {0, 3}

P:=seq(Or(x>p,x<p), p=discont(f,x));
                   P := Or(0 < x, x < 0), Or(3 < x, x < 3)

sols:=[solve( {g, P}, {x} )];
                      [                   /    3   1   (1/2)\ ]
              sols := [{3 < x}, {x = 1}, { x = - - - 13      }]
                      [                   \    2   2        / ]

eval(f, sols[3]);
                      /3   1   (1/2)\   /  3   1   (1/2)\
                    ln|- - - 13     | ln|- - - - 13     |
                      \2   2        /   \  2   2        /

evalf(simplify(%));
                                -11.29706355

It's also possible that options explicit and allsolutions might sometimes be useful/needed.

@Kitonum 

It's not very clear to me why solve finds this form more tractable. So it may or may not be of more general use.

kernelopts(version);

          Maple 2018.0, X86 64 LINUX, Mar 9 2018, Build ID 1298750

restart;
solve( evalc(Im(sqrt(x))), {x} );

                                  {0 <= x}

restart;
solve( evalc(Im(sqrt(x))), x );

                           RealRange(0, infinity)

Of course I realise that one may still run into the original bug reported here in the Question. I submitted a bug report against that earlier today (also with a shorter example of it).

My point is that even though using evalc may sometimes emit an error doesn't mean it couldn't be a principal line of attack.

But it may be difficult to delineate the trigger of similar buggy cases.

restart;
f:=sqrt(x)*ln(y):

g:=evalc(Im(f)): lprint(g);
(1/2)*abs(x)^(1/2)*(1-signum(x))*ln(abs(y))+(1/2)*abs(x)^(1/2)*(1+signum(x))*(1/2-(1/2)*signum(y))*Pi

solve(g, {x,y});
Error, (in unknown) invalid input: Check expects its 1st argument, sys, to be of type SubSystem, but received {-4*Pi, 0 <= x, 0 <= y, y <= 0}

h:=simplify(g) assuming real: lprint(h);
-(1/4)*x^(1/2)*(Pi*signum(y)-Pi+2*ln(abs(y)))*signum(x)^(1/2)-(1/4)*abs(x)^(1/2)*(Pi*signum(y)-Pi-2*ln(abs(y)))

solve(h, {x,y});
      {x = x, 0 < y}, {y = 1, x <= 0}, {x = 0, y < 0}, {y = -1, x < 0}

@ianmccr This is not a topic of advanced usage. Its just the basics of the syntax for calling a procedure (a.k.a. making a procedure call, or applying a procedure).

Your examples are not intrinsically connected with the subtlety of special evaluation rules.

Btw, a procedure which gets passed around (or returned, or applied to some arguments) without its being assigned to a name (ie. the procedure rather than the name is passed around, or used) is sometimes called an anonymous procedure.

nb. If one is trying to make a technical distinction then I don't much approve of using terms like "executing" a procedure, or "evaluating" a procedure, to denote applying a procedure.

 

First 248 249 250 251 252 253 254 Last Page 250 of 592