mmcdara

6149 Reputation

17 Badges

9 years, 72 days

MaplePrimes Activity


These are replies submitted by mmcdara

@Kitonum  @dnaviaux

Funny, we didn't understand the same thing at all.I
I didn't think of "display" in terms of "graphic display".

restart:
interface(version)
# one solution

N := 10:
V := Vector(combinat:-randperm([$1..10])):
Vector(N, [seq(V__||i = V[i], i=1..N)]);


# another one

seq(printf("%-5a = %2a\n", V__||i , V[i]), i=1..N)


# This one takes inspiration from Kitonum's reply

V := Vector([2, 1]):  
if V::Vector then V := convert(V, list) end if:

if numelems(V) = 2 then
  plots:-arrow(
     V,
     scaling=constrained, 
     caption=cat("first:", sprintf("%a", V[1]), "\nsecond:", sprintf("%a", V[2])),
     captionfont=[times,16]
  ):
elif numelems(V) = 3 then
  plots:-arrow(
     V,
     scaling=constrained, 
     caption=cat("first:", sprintf("%a", V[1]), "\nsecond:", sprintf("%a", V[2]), "\nthird:", sprintf("%a", V[3])),
     captionfont=[times,16]
  ):
else
  error "can't plot in spaces of dimension larger than 3, please call Dr Who to fix this"
end if;




Maple Worksheet - Error
Failed to load the worksheet /maplenet/convert/display.mw .

 

display.mw

 

@Carl Love 

In 

DEGREES(HW): ^2 ^1 ^0 ^1 ^0 

the first exponent (^2) is the total degree of the (first) monomial.
Next 

^1 ^0 ^1 ^0 

is just the power of each variables (according to the order in which thay are declared as NAME) in the first monomial.
Thus 

DEGREES(HW): ^2 ^1 ^0 ^1 ^0 

means "the first monomial equals a^1*b^0*x^1*y^0 = a*x"

It seems that this first exponent is redundant 
See for instance
 

dismantle(a^1*x^2*b^3*y^4+a^4*x^3*b^2*y^1);

POLY(6)
   EXPSEQ(5)
      NAME(4): a
      NAME(4): b
      NAME(4): x
      NAME(4): y
   DEGREES(HW): ^10 ^4 ^2 ^3 ^1 
   INTPOS(2): 1
   DEGREES(HW): ^10 ^1 ^3 ^2 ^4 
   INTPOS(2): 1

 

@nm 
 

Do not tell on user.wolfram.com that Mathematica has a feature Maple hasn't, or this will appear on :-)
compare-mathematica-and-maple.html   :-)

More seriously, your reply is very interesting for it could give the Maplesoft Development Team the idea of doing the same thing in Maple.

Thank you nm.

@acer 

"Certainly there is more work to be done in order to manipulate structures returned by ToInert": I agree.

For information:
"I don't know what kinds of expression you wish to cover":

my original issue was to represent graphically (and manipulate) an expression such that F0(X0, F) where X0 is a sequence (possibly void) of variables and F a sequence (possibly void) of functions of the form of F0.
A simple example  is F0(x, F1(y, z), F2(F3(u, v, w), F4(p, q)))  (note that all the set of variables in each Fn are mutually disjoint)
Here is the associated tree graph:

restart
with(GraphTheory):
G := Graph({{F__4, p}, {F__4, q}, {F__3, u}, {F__3, v}, {F__3, w}, {F__2, F__3}, {F__2, F__4}, {F__1, y}, {F__1, z}, {F__0, x}, {F__0, F__1}, {F__0, F__2}}):
DrawGraph(G, style=tree, root=F__0)


This kind of graph and/or expression appears in probability theory: for the example above, p, q, u, v, w, x, y, z are random variables and F0..F4 are functions named "copulas" that define some kind of dependency between the variables they operate on.
Copula_(probability_theory)
Such arborescent structure is called a "hierarchical copula" (HC for short, while hierarchical archimedian copula are far more common), see for instance here
HAC.pdf

I have developped a package to "sample from a HC" (at the image of the R package this last reference refers to) but I wanted two things:

  1. that the user can define directly (as in R) the dependency structure with a formula (the only thing that really matters to sample the HC) but taht he can visualize this formula graphicallyt (which is more easily understandable for other persons);
  2. or, that the user can define the structure graphically (I use Maplets to do this), which implies that  I have to construct the formula from the graph.
     

In fact I have to solve a direct (formula --> Graph) and a reverse (Graph --> formula) problem.
I just came back recently on the first point above and "extended" what I had done to other more general mathematical expressions.
That's why I spoke of this work as "an amusement or an exercice".

 

@acer 

Hi, 
Yes, I was aware of the existence of dismantle but I find that its results are less readable than a real graph... at least as long as the latter is not too big!
As you talk about dismantle, could you explain me this phrase (Maple 2015) taken from it's help page :
"...To produce an easily-manipulate inert form of a Maple structure, use ToInert instead."

When I read it first time, I had hoped that Maple functions did exist to manipulate the output of ToInert: it's probably because I didn't find any that I made this post. Nevertheless do not take it too seriously, it's more  fun or an exercise than anything else (although I think graphical visualization of a mathematical expression is pretty cool).

@nm 

This procedure covers a wider range of situations... maybe it will be enough for you?

F.mw

To overcome possible reading problems the code is here

F := proc(expr)
local nf, df, Expr:
nf := map(numer, [op(expr)]);
df := map(denom, [op(expr)]);
Expr := zip((u, v) -> if patmatch(u, y^b::rational*a::anything,'la') then 
                        [a, b] =~ eval(eval([a/v, b], la), y=__y)
                      else 
                        u 
                      end if, 
                      nf, df
        );
map(u -> if type(u, list) then
           if nops(rhs(u[2])) = 1 then
             eval(a*y^b, u)
           else 
             eval(a*y^b*__y^c, [u[1], b=iquo(op(rhs(u[2]))), c=irem(op(rhs(u[2])))/op(2, rhs(u[2]))])
           end if:
         else
           u
         end if,
         %
);
return eval( convert(add(%), horner, y), __y=y)
end proc:

expr := x+y^2*h(x, y)+y^3*f(y)*sin(x)+3*x*y^(7/2);
F(expr)
               2            3                    (7/2)
          x + y  h(x, y) + y  f(y) sin(x) + 3 x y     

            /          /                 (1/2)  \  \  2
        x + \h(x, y) + \f(y) sin(x) + 3 y      x/ y/ y 


expr := y^2+y^(8/3)/(y+1);
F(expr)
                                (8/3)
                           2   y     
                          y  + ------
                               y + 1 

                        /     (2/3)\   
                        |    y     |  2
                        |1 + ------| y 
                        \    y + 1 /   


 

@nm 

ok, thanks for the feed back

@acer 

Initially, I just wanted details on I , and as is often the case on Mapleprimes, the subject has slightly shifted to other considerations.
Nevertheless thank you for your comments.
 

@Carl Love 

Sorry, I've missed the "where a and are manifest real constants" point.
I understand now your what you have said.
 

@Carl Love 

Thank you for this veri documented reply.

Concerning Complex I had indeed observed it was a little peculiar while doing this:

ToInert(Complex(a, b)); #1
    _Inert_FUNCTION(_Inert_NAME("Complex", _Inert_ATTRIBUTE(

      _Inert_NAME("protected", 

      _Inert_ATTRIBUTE(_Inert_NAME("protected"))))), 

      _Inert_EXPSEQ(_Inert_NAME("a"), _Inert_NAME("b")))
ToInert(irem(a, b)); #2
      _Inert_FUNCTION(_Inert_ASSIGNEDNAME("iquo", "PROC", 

        _Inert_ATTRIBUTE(_Inert_NAME("protected", 

        _Inert_ATTRIBUTE(_Inert_NAME("protected"))))), 

        _Inert_EXPSEQ(_Inert_NAME("a"), _Inert_NAME("b")))
ToInert(f(a, b)); #3
      _Inert_FUNCTION(_Inert_NAME("f"), 

        _Inert_EXPSEQ(_Inert_NAME("a"), _Inert_NAME("b")))

While ToInert says Complex is a function, it's name is not assigned (like irem, cos, ln, ...) and looks like f in #3.
But at the same time Complex differs from f  because it has _Inert_ATTRIBUTE(s)  and protected _Inert_NAME(s), just as irem in #2.
It made me feel like Complex was someting in between an "usual" function (irem, cos, ln, ...) and a "user" one (f)

@nm 

Thank you nm.

@acer @vv 

So it is an alias!

Are there other "predefined" aliases in Maple?

 

Thank to both of you
To acer : I think my 2015 is too older for Help page for topic I, but I followed the link and it is useful.
 

@Carl Love 

"Am I so unlucky?"  
Not at all, it's an issue I have to deal with on a regular basis (Maple 2015, Mac OSX Mojave)
This occured for the first time about a year ago.
I remember Tom Leslie having the same problem suggested that forcing interface(rtablesize=10) at the top the of the worksheet could sometimes helped.

From my experience

  • this trick doesn't seem to be a panacea;
  • the problem is not related to the size of data the worksheet generates; for instance the file wich contains this simple lines cannot be displayed
    restart
    interface(version);
    I__D__moy := -(-1 + R)*I__L__moy;
    collect(expand(I__D__moy), I__L__moy):
    map(sort, %, R, ascending);
    Maple Worksheet - Error
    Failed to load the worksheet /maplenet/convert/arrange.mw .

    Download arrange.mw

     

  • it is more common when you use some packeges such that Statistics, plots/plottools


This of course doesn't help

@Kitonum 

Ok, I was fooled by "solve" and, out of overconfidence, didn't check his result.
Good point, I piteously vote up 

@Robert Matlock 

Could it be that you want to make sure that the numerical scheme used by dsolve is accurate enough for your purpose?
That is that you would want to insure the truncation error is small enoudh or that the modified equation (the continuous one the numerical scheme really solves) is close enough to the continuous equation you want to solve?

If it is so, there exist mathematical techniques to obtain either the truncation error or the modified equation.
Truncation error is the most classic and all math textbooks about numerical solving of ODEs contain (for toy problems, mainly linear) its expression for a lot of schemes (most of them avaliable in dsolve).

For more complicated case (ie. nonlinear odes) one often use the Method of Manufactured Solution (MMS).
The attached file gives an example.
Applying this method to the differential operator of your own problem, whilecarefully chosing several manufactured solutions, should help you to discriminate amid the "good" and "bad" numerical schemes and set the tolerances or grid size values.
 

restart:

ode := diff(u(x), x$2)*(1+diff(u(x),x)^2) + u(x)*diff(u(x),x) = 1/(1+15*x^2);
bv  := u(0)=1, u(1)=1/2;

sol := dsolve({ode, bv }, numeric);

plots:-odeplot(sol, [x, u(x)], x=0..1);

# Is this a reliable approximation of the true solution?

(diff(diff(u(x), x), x))*(1+(diff(u(x), x))^2)+u(x)*(diff(u(x), x)) = 1/(15*x^2+1)

 

u(0) = 1, u(1) = 1/2

 

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(8, {(1) = .0, (2) = .12826886497223997, (3) = .26113777569397906, (4) = .40154761998091554, (5) = .5494211741595527, (6) = .7027818747480806, (7) = .8595681924875884, (8) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(8, 2, {(1, 1) = 1.0, (1, 2) = -.8414172775411444, (2, 1) = .9007826428616932, (2, 2) = -.7082299965860446, (3, 1) = .8144476066224402, (3, 2) = -.5966120260026161, (4, 1) = .7371072611596353, (4, 2) = -.5095254313980015, (5, 1) = .6670249224566486, (5, 2) = -.4415965878326471, (6, 1) = .6035883533804827, (6, 2) = -.3880057661417157, (7, 1) = .546239605544721, (7, 2) = -.3452272743149798, (8, 1) = .5, (8, 2) = -.3143013506393865}, datatype = float[8], order = C_order); YP := Matrix(8, 2, {(1, 1) = -.8414172775411444, (1, 2) = 1.07812386883615, (2, 1) = -.7082299965860446, (2, 2) = .9589961107852607, (3, 1) = -.5966120260026161, (3, 2) = .7229274327553343, (4, 1) = -.5095254313980015, (4, 2) = .5303930579097911, (5, 1) = -.4415965878326471, (5, 2) = .3978675061925884, (6, 1) = -.3880057661417157, (6, 2) = .3069166956844302, (7, 1) = -.3452272743149798, (7, 2) = .2424437573985556, (8, 1) = -.3143013506393865, (8, 2) = .1999031726404807}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(8, {(1) = .0, (2) = .12826886497223997, (3) = .26113777569397906, (4) = .40154761998091554, (5) = .5494211741595527, (6) = .7027818747480806, (7) = .8595681924875884, (8) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(8, 2, {(1, 1) = .0, (1, 2) = 0.10890358515582652e-6, (2, 1) = -0.1273188731185599e-6, (2, 2) = 0.1433045594562052e-6, (3, 1) = -0.8927520188200564e-7, (3, 2) = 0.58257129736926296e-7, (4, 1) = -0.58336959823718046e-7, (4, 2) = 0.10428023193887122e-6, (5, 1) = -0.41584070894710373e-7, (5, 2) = 0.10452293110652422e-6, (6, 1) = -0.26811723670178546e-7, (6, 2) = 0.9748340266920589e-7, (7, 1) = -0.12413653369215185e-7, (7, 2) = 0.9099748924152847e-7, (8, 1) = .0, (8, 2) = 0.8633932517444991e-7}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[8] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.433045594562052e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 8, [u(x), diff(u(x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(8, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(8, 2, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[u(x), diff(u(x), x)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[8] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.433045594562052e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 8, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(8, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(8, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(0..0, {}), (3) = [x, u(x), diff(u(x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[u(x), diff(u(x), x)]'[i] = res[i+1], i = 1 .. 2)] catch: error  end try end proc

 

 

# MMS : Let's choose a function which verifies the bv

U__MMS := x -> 1-x/2

proc (x) options operator, arrow; 1-(1/2)*x end proc

(1)

# Obviously ode is not verified if one replace u by U__MMS.
# But if we replace its rhs by rhs__MMS computed this way

rhs__MMS := eval(lhs(ode), u(x)=U__MMS(x));

-1/2+(1/4)*x

(2)

# then the solution of

 ode__MMS := diff(u(x), x$2)*(1+diff(u(x),x)^2) + u(x)*diff(u(x),x) = rhs__MMS

(diff(diff(u(x), x), x))*(1+(diff(u(x), x))^2)+u(x)*(diff(u(x), x)) = -1/2+(1/4)*x

(3)

# should be U__MMS.
# Let's verify this

sol__MMS := dsolve({ode__MMS, bv }, numeric);

plots:-display(
  plots:-odeplot(sol__MMS, [x, u(x)], x=0..1, color=red),
  plot(U__MMS(x), x=0..1, style=point, symbol=circle, numpoints=20, color=blue),
  gridlines=true
);

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(8, {(1) = .0, (2) = .14285714285714282, (3) = .28571428571428564, (4) = .4285714285714285, (5) = .5714285714285715, (6) = .7142857142857144, (7) = .8571428571428572, (8) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(8, 2, {(1, 1) = 1.0, (1, 2) = -.5, (2, 1) = .9285714285714286, (2, 2) = -.5, (3, 1) = .8571428571428572, (3, 2) = -.5, (4, 1) = .7857142857142857, (4, 2) = -.5, (5, 1) = .7142857142857142, (5, 2) = -.5, (6, 1) = .6428571428571428, (6, 2) = -.5, (7, 1) = .5714285714285714, (7, 2) = -.5, (8, 1) = .5, (8, 2) = -.5}, datatype = float[8], order = C_order); YP := Matrix(8, 2, {(1, 1) = -.5, (1, 2) = .0, (2, 1) = -.5, (2, 2) = 0.5551115123125783e-17, (3, 1) = -.5, (3, 2) = 0.11102230246251566e-16, (4, 1) = -.5, (4, 2) = -0.22204460492503132e-16, (5, 1) = -.5, (5, 2) = -0.22204460492503132e-16, (6, 1) = -.5, (6, 2) = .0, (7, 1) = -.5, (7, 2) = .0, (8, 1) = -.5, (8, 2) = .0}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(8, {(1) = .0, (2) = .14285714285714282, (3) = .28571428571428564, (4) = .4285714285714285, (5) = .5714285714285715, (6) = .7142857142857144, (7) = .8571428571428572, (8) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(8, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0, (8, 1) = .0, (8, 2) = .0}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[8] elif outpoint = "order" then return 2 elif outpoint = "error" then return HFloat(-0.0) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 8, [u(x), diff(u(x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(8, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(8, 2, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[u(x), diff(u(x), x)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[8] elif outpoint = "order" then return 2 elif outpoint = "error" then return HFloat(-0.0) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 8, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(8, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(8, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(0..0, {}), (3) = [x, u(x), diff(u(x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[u(x), diff(u(x), x)]'[i] = res[i+1], i = 1 .. 2)] catch: error  end try end proc

 

 

# error

plots:-odeplot(
  sol__MMS, [x, u(x)-U__MMS(x)], x=0..1,
  color=red, gridlines=true,
  labels=["x", "error"], labeldirections=[default, "vertical"]
)

 

 


 

Download MMS.mw

 

 

First 65 66 67 68 69 70 71 Last Page 67 of 125