Carl Love

Carl Love

28025 Reputation

25 Badges

12 years, 313 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

If I define w(x) = int(y*u(y)^3, y= 0..x) and w(1) = J, then your IVP becomes the BVP system

eq1:= diff(w(x),x) = x*u(x)^3, w(0)=0, w(1)=J;
eq2:= diff(u(x),x) = exp(x) + (2*exp(3)-1)/9 + J, u(0)=1; 

Hopefully the derivation of that is obvious; if not, let me know. The system has a parameter, J. Maple will numerically solve BVPs with parameters provided that there's an additional BC for each parameter. As you can see, we have 3 BCs and the differential order is 2; so, we're good to go.

dsn:= dsolve({eq1, eq2}, numeric);
Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

Uh oh, we hit a roadblock. This is a very common error message from Maple's BVP solver. When you get this error (note the word initial), the first thing to try is a kind of bootstrapping called a continuation parameter. This is an additional parameter C that varies continuously from 0 to 1 that you place anywhere in the system (usually as a coefficient) such that when C=1, you get the original BVP, and when C=0, you have a relatively simpler BVP (see ?dsolve,numeric,BVP). It may take some guessing and a few trials to find a place to put C that works. (You can also put C in multiple places at once.)

In this case, what works is changing the BC w(1)=J to w(1)=J*C:

eq1:= diff(w(x),x) = x*u(x)^3, w(0)=0, w(1)=J*C;
eq2:= diff(u(x),x) = exp(x) + (2*exp(3)-1)/9 + J, u(0)=1;
dsn:= dsolve({eq1, eq2}, numeric, continuation= C);
plots:-odeplot(dsn, [x,u(x)], x= 0..1); 

It works, and the resulting plot and a deeper analysis of the numeric values shows complete agreement with Mariusz's symbolic solution.

If you need to continue the solution for u(x) beyond the interval 0..1, that can be done with a little extra work.

Okay, here is a procedure to do it:

VuReduce:= proc(e)
uses P= Physics;
local old:= P:-Expand(e), r;
   do
      r:= (P:-Expand@evalindets)(
         old,
         specfunc(P:-`*`),
         proc(`*`)
         local p,p1,U,V;
            if
               membertype({v[nonnegint], 'P:-`^`'(v[nonnegint], posint)}, `*`, 'p')
               and nops(`*`) > p
               and (U:= op((p1:= p+1), `*`))
                  ::{u[nonnegint], 'P:-`^`'(u[nonnegint], posint)}
            then
               V:= op(p,`*`);
               subsop(
                  p= Vu(
                     `if`(V::v[nonnegint], op(V), op([1,1],V)), # "a"
                     `if`(U::u[nonnegint], op(U), op([1,1],U)), # "b"
                     `if`(V::v[nonnegint], 1, op(2,V))          # "c"
                  ),
                  p1= `if`(U::u[nonnegint], NULL, P:-`^`(op(1,U), op(2,U)-1)),
                  `*`
               )
            else
               `*`
            fi
         end proc
      );
      if r = old then return r fi;
      old:= r
   od
end proc:

Call it via

VuReduce(e)

where e is any expression (or container of expressions) that you want to apply the transform to. Please test this. I made it so that it will  keep applying the transform as long as the expression keeps changing under it. Hopefully this won't lead to an infinite loop, but I can't know that without knowing the details of Vu. If it's a problem, I can easily include code to detect cycling. Note that I apply Physics:-Expand to both the original and transformed expressions.

If you have any trouble with this, please post a worksheet showing the erroneous output or error message and an explanation of what you think the result should've been.

I don't know if you're in the habit of using either 1D or 2D Input. I only use 1D. I haven't tested the above procedure entered as 2D Input. If this makes any difference, it would only affect the entry of the procedure itself. The input form of the expressions or the procedure's call won't matter. I'm not expecting there to be any differences, but there are many unfortunate surprises with how 2D parses the code.

I've never worked with the Physics package before, so there may be some easier way to do this that I missed. Nonetheless, the above should work with reasonable efficiency.

You could do it like this:

Trunc:= proc(F::{procedure, `+`}, odr::posint:= 2, v::list(name):= [x,y,z])
local f:= `if`(F::procedure, F(v[]), F);
   if not f::`+` then
      error "Can't truncate 1-term expression"
   else
      select(q-> degree(q,v) <= odr, f)
   fi
end proc:

 

As hinted at by Tom, the only purpose of ArrayTools:-Alias is to create a reshaped, re-indexed, resized, and/or re-ordered secondary reference to an rtable (an Array, Vector, or Matrix) or to some subpart of it. If you apply Alias directly to a Matrix(...) command, then there's no primary reference to that Matrix, so there can be no point to creating a secondary reference (although it will allow you to do so). A primary reference would be created by first assigning the result of Matrix(...) to a variable.

 

You have not declared Wn to be an Array. When an assignment is made to an indexed undeclared variable, a table (but not an rtable) is created. But table indexing works differently than Array indexing in that index ranges such as 0..1 are not respected. In your case, Wn ends up with only 3 indices---[0..1, 0..1, 1], [0..1, 0..1, 2], and [0..1, 0..1, 3]---so the ranges are literally part of each individual index. You can inspect the end result of your loop with eval(Wn).

Representing  (1 - 1e-18) requires 18 significant digits, whereas representing 1e-18 requires only 1 significant digit. So I find the following way much easier than all the others presented in this thread, requiring no increase in Digits from its default. I apply the following obvious identity:

Quantile(Normal(mu,sigma), 1-x) = mu - sigma*Quantile(Normal(0,1), x).

Thus

Quantile(Normal(0,1), 1 - 1e-18) = -Statistics:-Quantile(Normal(0,1), 1e-18);
      8.75729034878321

The following procedure will apply the rule only to those lns that occur as a result of the integration:

REALINT:= proc(f::algebraic)
local OrigLn:= indets(f, specfunc(ln));
   evalindets(
      int(args), 
      specfunc(ln), 
      L-> `if`(L in OrigLn, L, ln(abs(op(L))))
   )
end proc:

Example:
REALINT(1/(x*ln(x)), x);

If you're willing to have Maple tacitly make the necessary assumptions that'll get you that most simplification, you can use simplify(..., symbolic). Although I've never seen this occur in a practical situation, it is possible to construct expressions where the necessary assumptions are contradictory, and thus the simplification is not valid for any possible valuations.  See ?simplify,symbolic.

Manipulating square roots is difficult. Here's one tip: Almost all sqrt(...) are converted to (...)^(1/2), so searching expressions for sqrt won't work.

Here's a procedure to integrate recursively until either there's a term with no derivative or we get stuck with an unevaluated integral:

MyOdeInt:= proc(ODE::algebraic, Z::name:= NULL)
local 
   ode:= frontend(expand, [ODE]), #Distribute, but don't expand functions.
   #Find the integration variable:
   z:= `if`(Z=NULL, map2(op, -1, indets(ode, 'specfunc'(diff)))[], Z),
   r 
;
    if ode::`+` then #There's more than 1 term.
       if remove(hasfun, ode, diff) <> 0 then return ode fi 
    #single-term case:
    elif not hasfun(ode, diff) then return ode 
    fi;
    r:= int(ode, z);
    `if`(hasfun(r, int), r, thisproc(r))
end proc:

And you'd use it:

MyOdeInt(odetemp);

I tested on higher orders:

MyOdeInt(diff(odetemp, z$9));

which returned the same thing as the first command.

At the beginning of your worksheet, do

with(Units); with(Standard);
UseUnit('J/(mol*K)'):
UseUnit('J/(kg*K)'):

Then computations whose results can be expressed in those units will be. Computations include arithmetic, evalf, and combine(..., units). The raw extraction of R from ScientificConstants will stay the same, but you can force it to the preferred units with simply combine(..., units) rather than needing convert.

HFloat stands for "hardware float". These are the 64-bit representations of floating-point (or decimal) numbers that are nearly universally used by computer chips. (They're called "double precision" in many languages.) Included are representations of infinity, -infinity, and undefined.

For most purposes, you can think of HFloat(infinity) as simply infinity. If it appears as a result of integration, it indicates some type of divergence issue. Perhaps the integral is truly divergent. If it appears multiplied by some coefficient expression C, and C could possibly be 0, then you may have an indeterminate form, which'll require deeper analysis. There's no simple way to deal with this. It requires case-by-case analysis.

The results from solve should be put into a set. Then, iterate over the members (my s) of the set without regard for the index numbers (your k). Like this:

f:= x^2+2*x*y+2*y^3;
V:= [x,y];
gr:= VectorCalculus:-Gradient(f,V);
soln:= {solve({seq(gr)}, {V[]})};
(H,HD):= VectorCalculus:-Hessian(f, V, 'determinant');
seq(eval([H,HD], s), s= soln); #Change "s=" to "s in" if you like. 
#or
map2(eval, [H,HD], soln);

 

The reason for your error is that nops can't be applied to a sequence (unless it's a trivial sequence with 1 member)[1]:

S:= 3,5,7,9: #a sequence
nops(S); #error
nops([S]); #make S a list

      4

So, it'd be a rare situation that nops(op(X)) would work at all! It could only work if nops(X) = 1. Also, note that nops(X) is by definition (see ?nops) the number of entries of [op(X)].

[1]: Indeed, there's only a handful of cases where a nontrivial sequence can be used as a single argument to a procedure, and those are all built-in procedures.

It seems to be a bug. Here's a workaround:

J1:= int(exp(-t)/(1-t), t= 0..1-epsilon) assuming epsilon > 0, epsilon < 1;
J2:= int(exp(-t)/(1-t), t= 1+epsilon..infinity) assuming epsilon > 0;
simplify(J1+J2);
limit(simplify(J1+J2), epsilon= 0, right);

Yielding Ei(1)/exp(1), a real value. If the limit is incorrectly taken from the left, a related nonreal value is returned. Perhaps this has something to do with the bug.

Here's an implementation of procedure evalpw, as mentioned by Joe. I believe that I've covered every possibility; if not, let me know.

evalpw:= proc(PW, V::{`=`, {set,list}(`=`)})
description 
   "reduces a piecewise expression by considering its conditions" 
   " under evaluation without evaluating the branches"
;
option `Author: Carl J Love <carl.j.love@gmail.com> 3-Aug-2018`;
(* 
The 1st arg can be anything. If it's not a piecewise,
then it's returned unchanged. The 2nd arg is exactly like 
eval's 2nd arg, but it's only used to temporarily evaluate 
the conditions.  

The conditions are evaluated *in order* using PiecewiseTools:-Is
(undocumented, but its code is easy to read). If any evaluates
true, its branch is immediately returned. Any that evaluate false are 
removed from further consideration along with their branches. If all
evaluate false, then a default or "otherwise" clause is returned.
Otherwise, a reduced piecewise expression is returned.
*)
local  
   pw:= [op(PW)], n:= nops(pw), 
   k, #index to operands of PW  
   b, #boolean evaluation of condition
   R:= table() #undecidable conditions and their branches
;
   if not PW::'specfunc'(piecewise) then return PW fi;

   #Iterate over the conditions:
   for k by 2 to n-1 do
      b:= PiecewiseTools:-Is(eval(pw[k], V));
      if b = true then return pw[k+1] fi;
      if b <> false then (R[k],R[k+1]):= (b,pw[k+1]) fi;
   od;

   if numelems(R) = 0 then #No true, all false
      #If n::odd, use "otherwise", else use default.
      #Note that the default may be an index on 'piecewise' 
      #(it's an undocumented feature).
      return `if`(n::'odd', pw[n], op(0,PW)(false, ``))
   fi; 
   
   #Build residual piecewise:
   if n::'odd' then R[n]:= pw[n] fi;
   op(0,PW)(entries(R, 'indexorder', 'nolist'))         
end proc:

 

Make the first line of code of your worksheet:

restart;

This line should be in its own execution group.

If that doesn't solve your problem, then you'll need to attach your worksheet to a post, which you should usually do anyway when asking a Question. Use the green up-arrow on the editor's toolbar.

First 166 167 168 169 170 171 172 Last Page 168 of 395