Carl Love

Carl Love

27626 Reputation

25 Badges

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

MaplePrimes Activity


These are replies submitted by Carl Love

@Thomas Richard The help page (?Student,NumericalAnalysis,MatrixDecomposition) phrase "... using the Crout factorization algorithm for tridiagonal matrices" is poorly worded and hence a bit misleading. The Crout algorithm itself is not restricted to tridiagonal matrices; rather, the help page is trying to say that this command's implementation of that algorithm is limited to tridiagonal matrices. See the subsection "LU Crout decomposition" of the Wikipedia page "LU decomposition".

With a tiny additional step, you can convert your LDU decomposition to a Crout: Set L:= L.D. However, by using transposes as shown in my Answer, the slight extra computation of factoring the from can be avoided.

@nm For some reason that I don't understand, the assumption that the argument of ln is positive isn't enough. If you also assume _C1 < 0, it'll work.

odetest will also work for an explicit separable solution by using the right assumptions. I'm just pointing this out because it's true, not because I think the separable solution is better. I understand why the d'Alembert solution is better.

restart:

ode:= y(x) = x + 3*ln(diff(y(x), x)):

It's not necessary to use an initial condition for this. I just do it to ensure that the constant of integration

is the same in my two solutions. However, note that my initial condition is completely arbitrary.

ic:= y(x0)=y0:

ivp:= {ode, ic}:

Solution by dsolve(..., [separable]) and verification with odetest:

sol1:= dsolve(ivp, [separable]);

y(x) = -3*ln(-exp(-(1/3)*x0)+exp(-(1/3)*y0)+exp(-(1/3)*x))

assmp:= sol-> ((x,x0,y0)::~real, op(indets(sol, specfunc(ln))[]) >~ 0):

odetest(sol1, ivp) assuming assmp(sol1);

{0}

"By-hand" solution via Calc II separable technique:

eq1:= expand(eval(solve(ode, diff(y(x),x)), y(x)= y));

exp((1/3)*y)*exp(-(1/3)*x)

(ys,xs):= selectremove(has, eq1, y);
solgen:= solve(int(1/ys, y) = int(xs, x) + C, y);

exp((1/3)*y), exp(-(1/3)*x)

3*ln(-3/(-3*exp(-(1/3)*x)+C))

sol2:= (simplify@eval)(
    y(x) = solgen,
    solve({rhs(ic) = eval(solgen, x= op(lhs(ic)))}, C)
);

y(x) = 3*ln(1/(-exp(-(1/3)*x0)+exp(-(1/3)*y0)+exp(-(1/3)*x)))

odetest(sol2, ivp) assuming assmp(sol2);

{0}

Verify equivalence of solutions:

simplify(rhs(sol1-sol2)) assuming assmp(sol2), assmp(sol1);

0

 

Download DsolveSeparable.mw

@Andiguys Yes, you can include any number of inequalities in the first argument to solve. Using them has no effect on the number-of-equations-versus-number-of-variables issue. (You can't use inequalities with eliminate.) For example,

sol:= solve({tau = `&tau;opt`, lambda = `&lambda;opt`, Rem > Rom}, {w, Clm});

(`&tau;opt` is the correct way to write tauopt in 1D input (Maple Notation); I incorrectly put tauopt in the Answer above. This comment doesn't matter if you're using 2D Input, which seems to be the case.)

This type of solution (see ?solve,parametric and ?SolveTools,Parametric) very often takes a long time. The solution is returned as a piecewise structure that is often too wide to be easily read. In that case, it helps to use lprint:

lprint(sol);

The returned form will be 

piecewise(condtion[1], [solution[1]], ..., condition[n], [solution[n]], [otherwise_solution]),

where each condition is an inequality, an equation, or several such in disjunctive normal form, i.e., 

Or(And(ineq[1][1], ..., ineq[1][k[1]]), ..., And(ineq[n][1], ..., ineq[n][k[n]])),

and otherwise_solution is the solution (often empty) if all the conditions are false. Note that each solution is presented as a list, even if it's a singleton. Thus, the kth condition can be extracted via op(2*k-1, sol) and its corresponding solution via op(2*k, sol)[]. The otherwise_solution can be extracted via op(-1, sol)[].

There is an asymmetry between your procedures and above that I find odd. I'd like you to check whether it's correct. In B, the term b2*(s2+h) occurs in in both factors of the denominator, and the similar term b1*(s1+v) also occurs once. In A, there are three distinct such terms: b1*(s1+h)b2*(s2+v), and b1*(s1+v). So, are those correct?

I confirm that I frequently encounter this bug in all Maple versions starting with Maple 2021, using either Windows 10 or Windows 11. This is not a rare bug for me. It happens roughly 30% - 50% of the time that I try to use a context menu.

@C_R The use of ? like an alphabetic letter in symbols is described in the very first sentence of help page ?nameI tend to use it as the last letter of the name of a procedure (including type checks) that returns true or false.

@zenterix You are correct: The 2nd argument to simplify will only be interpreted as side relations if it is a set or a list.

 

@C_R My guess is that FileTools hasn't kept up with file system changes made for Windows 11. The problem is that there are secret system "files" that are listed in directories, yet they don't "exist" by whatever criterion is used to determine existence. I don't know if this is new to Windows 11.

Anyway, the problem is easy to work around: We only check files and directories for which both FileTools:-Exists and FileTools:-IsReadable return true. Here's some code for it:

AllFiles:= module()
uses FT= FileTools;
local
    Q, dirs, DirsRead, files, Selected, Files, CurDir, time0, time1,   
    Readable?:= f-> FT:-Exists(f) and FT:-IsReadable(f), 
    Report:= proc(RepIntv)
        Files+= nops~([files, Selected]);
        if (time1:= time[real]()) - time0 > RepIntv then
            time0:= time1;
            printf(
                "Current dir: %s\n"
                "Dirs checked: %-7d  Dirs stacked: %-7d  "
                "Files checked: %-7d  Files selected: %-7d\n\n",
                CurDir, DirsRead, numelems(Q), Files[]
            )
        fi
    end proc,
    ModuleApply:= proc(base::string:= "", ext::string:= "mw", {RepIntv::nonnegint:= 15})
    local Wanted?:= f-> FT:-Extension(f) = ext;
        Q:= DEQueue(base);
        Files:= [0,0];
        time0:= time[real]();
        [
            for DirsRead do
                (dirs, files):= selectremove(
                    FT:-IsDirectory, 
                    select(Readable?, FT:-ListDirectory((CurDir:= back(Q)), 'absolute'))
                );
                Selected:= select(Wanted?, files);
                pop_back(Q);  Q,= dirs[];
                if RepIntv <> 0 then Report(RepIntv) fi;              
                Selected[]
            until empty(Q)
        ]
    end proc
;
end module
:

F:= AllFiles():

Current dir: C:\\Windows\WinSxS\Temp\InFlight\50b12bd18ca6da016e0700007868d418\wow64_microsoft-windows-s..ity-netlogon-netapi_31bf3856ad364e35_10.0.22621.3085_none_bda3256405ca58e6
Dirs checked: 23640    Dirs stacked: 12102    Files checked: 11452    Files selected: 0      

Current dir: C:\\Windows\WinSxS\Manifests
Dirs checked: 50587    Dirs stacked: 10235    Files checked: 31094    Files selected: 0      

... many lines omitted ...

It takes about an hour for it to work through the 100,000+ "existing" and "readable" directories on my computer.

Can someone please check whether this code works on Linux and/or Apple?

All of the actual work can be done by the following code; everything else above is just to display the status updates:

AllFiles:= proc(base::string:= "", ext::string:= "mw")
uses FT= FileTools;
local
    Q:= DEQueue(base), dirs, files, 
    Readable?:= f-> FT:-Exists(f) and FT:-IsReadable(f),
    Wanted?:= f-> FT:-Extension(f) = ext
;
    [
        do
            (dirs, files):= selectremove(
                FT:-IsDirectory, 
                select(Readable?, FT:-ListDirectory(pop_back(Q), 'absolute'))
             );
             Q,= dirs[];
             select(Wanted?, files)[]
         until empty(Q)
    ]
end proc
:

 

I think that it would be easy to automate this by using FileTools:-ListDirectory to find all "*.mw" files on one's computer and then using eBookTools:-CreatePDF to make PDFs from them.

@salim-barzani When we say "Upload your worksheet", worksheet refers to the file that contains a record of what you see on your screen during your Maple session.

I notice that the article's solutions have no constant of integration. Did the article specify an initial condition? I think that the article's solutions are equivalent to the 3rd and 4th Maple solutions with _C1 set to 0, but it'll take some effort to verify that. I doubt that you'll ever be able to get Maple to arrange the solutions in piecewise form on its own, but you should be able to verify the correctness of the given piecewise.

The 1st and 2nd solutions that Maple gives are singular solutionsNote that they don't depend on xi, so they make diff(z(xi), xi) identically equal to 0.

@Andiguys A strict inequality is one of the form a < b or a > b; a non-strict inequality is one of the form a <= b or a >= b. Strict inequalities are not allowed as constraints for real-valued variables by numerical optimization software. This is essentially the same as saying that the feasible set must be a closed (but not necessarily bounded) set (in the sense of real analysis); however, I think that the error message comes from a purely syntactic (rather than mathematical) checking of the constraints.

Here's a simple thought experiment to show why strict inequalities aren't allowed:

  • What is the minimum value of x^2 given that x > 0?

You have an unspecified parameter F in the line defining Dp.

Your z seems to be an independent variable. You can only make a 2-D plot such as you show for a specific value of z, which you haven't given.

Your sample plot says M^2 = 2, 5, 7. Did you want or M^2 set to those values? either is easy.

I corrected about 2 dozen missing-multiplication-sign errors in your code, and made a few other stylistic changes, obtaining this:

h:= z-> piecewise(
   z <= d+1,   1,
   z <= d+4,   1-delta/2*(1 + cos(2*Pi*(z - 1 - 1/2))),
   z <= d+6,   1 
):
w0:= -c*h(z)^2/4 + 3/64*(b*c-4*a)*h(z)^4 + 19/2304*b*(b-4*a)*h(z)^6:
w1:= c/4 + (4*a-b*c)*h(z)^2/16:
w2:= (4*(b*c-4*a)-b*h(z)^2)/256:
w3:= b*(b-4*a)/2304:
a:= x4*S*Gr*sin(alpha)/(4*x1*x5):
b:= 1/Da + x3*M/(x1*(1+m^2)):
c:= Dp/x1:
Dp:= 96*x1/((6-b*h(z)^2)*h(z)^4)*(F + a*h(z)/24 - 11/6144*b*(b-4*a)*h(z)^8):
x1:= ((1-phi1)*(1-phi2))^(-5/2):
x2:= (1-phi2)*(1 - phi1 + phi1*Rs1/Rf) + phi2*Rs2/Rf:
x3:= shnf/sf:
x4:= (1-phi2)*(1-phi1+phi1*RBs1/RBf+phi2*RBs2/RBf):
x5:= khnf/kf:
shnf:= sbf*(ss2+2*sbf-2*phi2*(sbf-ss2))/(ss2+2*sbf+phi2*(sbf-ss2)):
sbf:= sf*(ss1+2*sf-2*phi1*(sf-ss1))/(ss1+2*sf+phi1*(sf-ss1)):
khnf:= kbf*(ks2+2*kbf-2*phi2*(kbf-ks2))/(ks2+2*kbf+phi2*(kbf-ks2)):
kbf:= kf*(ks1+2*kf-2*phi1*(kf-ks1))/(ks1+2*kf+phi1*(kf-ks1)):
params:= {
    RBs1= 8933*16.7e6,
    RBf= 1063*1.8e6,
    RBs2= 6320*18e6,
    kf= 0.492,
    sbf= 6.67e-1, ss2= 2.7e-8,
    sf= 6.67e-1, ss1= 59.6e6,
    ks2= 76.5, kf= 0.492, ks1= 401,
    phi1= 0.01, phi2= 0.02, alpha= Pi/4, m= 0.5, Da= 0.1, Gr= 5, 
    delta= 1, S= 0.5,  d= 1
}:
                                      
W1:= eval(w0+w1*r^2+w2*r^4+w3*r^6, params):

 

@nm If you prefer the prettyprinted style of subscripted constants, there's no reason to stop using them. But they are aliases, so you mustn't treat them like variables. In code, you should use _C1, etc.

What I said before about aliases in procedures actually applies to aliases in any code: The alias must be defined before the code that uses the alias is entered. I only explained this with respect to procedures earlier because that's the most likely scenario for making this mistake.

In your example that doesn't work, you do a restart and then enter some code that contains a name which will become an alias. But it can't already be an alias at the time that the code is entered because it's immediately after a restart. The fact that c__1 will become an alias during the execution of the code is irrelevant. All that matters is whether it's already an alias when the code is entered (formally speaking, when the code is parsed).

I concur: Use a .zip file.

First 9 10 11 12 13 14 15 Last Page 11 of 704