Preben Alsholm

13471 Reputation

22 Badges

20 years, 218 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

The year in the copyright will say something, I suppose:

op(3,eval(applyrule));
    `Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2005`

On the other hand I'm sure simplify existed before this copyright year:
op(3,eval(simplify));
     system, `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved.`

Take a look at
u:=unapply(U(0.15036/6.022e23,0.035,200,T),T);
# and examine
evalf[15](u~(X));
# You will notice small imaginary parts and smaller if 15 is replaced by e.g. 20.
Thus we conclude that those are due to roundoff.
A simple way to get rid of them is to take the real part as in
(Re@u)~(X);
## Thus you can do:
Statistics:-NonlinearFit(Re(U(0.15036/6.022e23,d,Theta,T)),X,Y,T,initialvalues=[Theta=200, d=0.035]);

The design of assume (and Physics:-Assume) is not made to work as eval when assumptions are equalities as in your example.
It is rather interesting though that 0 does emerge as the output from the first of these:

restart;
assume(a+b=0);
for f in [x->sqrt(x^2),sin,exp,ln,x->1/x] do simplify(f(a+b)) end do;
####
restart;
Physics:-Assume(a+b=0);
for f in [x->sqrt(x^2),sin,exp,ln,x->1/x] do simplify(f(a+b)) end do;
#### and also from the assuming version:
restart;
for f in [x->sqrt(x^2),sin,exp,ln,x->1/x] do simplify(f(a+b)) assuming a+b=0 end do;
#######
If you remove simplify in the 3 examples above the output doesn't change.
################# Use simplify/side relations instead:
restart;
for f in [x->sqrt(x^2),sin,exp,ln,x->1/x] do simplify(f(a+b),{a+b=0}) end do;
#Output 0 0 1,error as expected from ln(0).
## You can catch the errors:
restart;
for f in [x->sqrt(x^2),sin,exp,ln,x->1/x] do
  try
    res:=simplify(f(a+b),{a+b=0})
  catch:
    res:=undefined;
  end try;
  res
end do;

The description in the help page:
?dsolve,numeric,bvp

makes it clear that the methods used are finite difference schemes: Notice e.g. the options initmesh and maxmesh.
Starting with a mesh over the interval and an initial approximation to the discretized problem. That initial approximation can be given in the optional argument approxsoln.

The method used is not rkf45 as that is purely for initial value problems.

You should begin by replacing all square brackets i.e. [ ], by parentheses, ( ), since square brackets are used for other things e.g. for making lists:  [a,b,c] is a list and so is [a] .

Secondly I see missing multiplication signs (*) in e.g. tau_mafw, tau_mafo, tau_mafg.

You cannot use a name like S_wfi both as a function name as it is in S_wfi(x,t) and also as a constant as it is in e.g. IBC.
Maybe the appearance of S_wfi(x,t) is an error and should have been S_wf(x,t) ?

After those corrections (I included the replacement of S_wfi(x,t) with S_wf(x,t) ) the variables can be found by these 3 steps

indets({Eq1, Eq10, Eq11, Eq12, Eq2, Eq3, Eq4, Eq5, Eq6, Eq7, Eq8, Eq9},specfunc(diff));
op~(1,%);
vars:=% minus indets(%,specfunc(diff));

But you will have no success with
pdsolve({Eq1, Eq10, Eq11, Eq12, Eq2, Eq3, Eq4, Eq5, Eq6, Eq7, Eq8, Eq9, IBC},vars);

which doesn't surprise me at all.
If you leave out IBC then pdsolve will work for a long time, but I interrupted the computation since I didn't expect an answer anyway.



Down to basics I would do:

restart;
v1:=<1,0,0>: v2:=<0,1,0>: v3:=<0,0,1>:
u1:=<0,1,0>: u2:=<0,0,1>:
A:=<v1|v2|v3>;
#We need to find out if both of the equations A.x = u1 and A.x = u2 have solutions.
# That is the same as asking if  A.X = <u1|u2> has a solution X.
# Thus we form the augmented matrix:
T:=<A|u1|u2>;
## and then simply:
LinearAlgebra:-Rank(T) - LinearAlgebra:-Rank(A); #Yes if 0, No if greater than zero
## You could solve, but that is not necessary here.
#LinearAlgebra:-LinearSolve(A,<u1|u2>);

I have no problem in Windows 10 and with the 64 bit version of Maple 2016.1.

In Maple 2016 the two prints in each procedure produce the same, thus you get
                              0, 0
                              1, 1
                              0, 0
                              1, 1
In Maple 8 it is as you report for Maple V Release 4:

                                 0, 0
                                 0, 1
                                 0, 0
                                 1, 1
So it seems that assign treats environtment variables different than other variables in early versions of Maple.
The help page for Environment Variables in Maple V Release 4 (and in Maple 8) starts by saying:
     "Environment variables can be used in a simple assignment."

But the help page for Maple 2016 says the exact same thing.
In Maple V Release 4 and in Maple 8 the global version
_EnvX:=0; assign('_EnvX=1');
doesn't change _EnvX to 1 either as does Maple 2016.

In Maple 2016 all four produce the same output:

map((a::uneval,b)->'args',[a,b,c,d],1..4,x);
map((a,b::uneval)->'args',[a,b,c,d],1..4,x);           
map((a,b)->'args',[a,b,c,d],1..4,x);
map((a,b)->args,[a,b,c,d],1..4,x);

## The output is
      [a, 1 .. 4, x, b, 1 .. 4, x, c, 1 .. 4, x, d, 1 .. 4, x]

When you say the product is Maple V, which release is it?

The output is the same in Maple 12 as in Maple 2016, but Maple 8 returns ['a, 1 .. 4, x', 'b, 1 .. 4, x', 'c, 1 .. 4, x', 'd, 1 .. 4, x'] in the first case, and [a, 1 .. 4, x, b, 1 .. 4, x, c, 1 .. 4, x, d, 1 .. 4, x] in the other three.

The function ux has complex values, so plot can only plot the real and/or the imaginary parts:

plot([Re@ux,Im@ux], 0..10, labels=[x,"Re(u),Im(u)"],color=[green,red],caption="t=0.8");

I suppose that the Digits setting at 5 is a misprint and ought to have been at least 50 as your data appears to have that many digits.
At Digits:=50 I don't get anything like the graph you show, not for Re(u) nor for Im(u).
Raising Digits to 60 still doesn't produce that graph. There are huge values at x > about 7.
Raising Digits to 100 produced this:



The numerical inversion of the laplace transform is notoriously difficult, I understand.

nb/lhs(hh)*rhs(hh);

I never use patmatch, so I cannot be of much help there. Notice though that you don't have a square root at all.

indets(f1,radical); #A fourth root
nops(f1); # 3
whattype(f1); # `*` a product
op(f1); #The 3 operands in the product
op(2,f1); # We can extract the radical and the polynomial part from this factor:
rad,pol:=selectremove(type,op(2,f1),radical);
rad;
pol;
## You shouldn't forget the two other factors, though:
op(1,f1);
op(3,f1);

There seems to be a solve problem as well as a LinearSolve problem. The reason is not clear to me right now.
(See note 3 at the bottom for another simple way out).
But your system SYS2 with the boundary conditions bcs2 and bcs22 can be solved like this.

Define SYS2, bcs2, and bcs22 as you did.
Then do:

Hdiff:={diff(h1(theta),theta$4),diff(h2(theta),theta$3),diff(h3(theta),theta$4)};
#solve(SYS2,Hdiff); #Does it finish?
A,b:=LinearAlgebra:-GenerateMatrix(convert(SYS2,list),[op(Hdiff)]);
# Try
LinearAlgebra:-Determinant(A); # Pretty small
LinearAlgebra:-ConditionNumber(A); #Pretty large
#LinearAlgebra:-LinearSolve(A,b); #Does it finish?
## That could be the reason for solve and fsolve not to return in a reasonable time; see note at bottom, though.
## Instead of using LinearSolve, we just do
sol:=A^(-1).b;
sys:=[op(Hdiff)]=~convert(sol,list):
res:=dsolve({op(sys)} union bcs2 union {bcs22},numeric);
plots:-display(Array([seq(plots:-odeplot(res,[theta,h(theta)]),h=[h1,h2,h3])]));
###############
Note added. The problem for LinearSolve seems to have something to do with the vector b rather than with the matrix A:
This causes absolutely no problem:
LinearAlgebra:-LinearSolve(A,<x,y,z>);
sol:=eval(%,{x=b[1],y=b[2],z=b[3]});
#Then continue as before.
Note 2 added: I tried some freezing when using solve with not that much success:
Hsbs:=Hdiff=~freeze~(Hdiff);
Fus:=indets(SYS2,function) minus Hdiff;
Fsbs:=Fus=~freeze~(Fus);
SYS2F:=subs(Hsbs,Fsbs,SYS2);
sol:=solve(SYS2F,rhs~(Hsbs)): #Use colon since it is huge.
sys:=thaw(sol);
dsolve(sys union bcs2 union {bcs22},numeric); #Complains about singularity at left endpoint
dsolve(sys union bcs2 union {bcs22},numeric,method=bvp[midrich]); #Complains about singularity
## I wouldn't trust these last two "results".
###############################
Note 3: Freezing expressions of type `^`  helps greatly:
SYS2FF:=evalindets(SYS2,`^`,freeze);
sys:=thaw(solve(SYS2FF,Hdiff));# FAST
# Now proceed with
res:=dsolve(sys union bcs2 union {bcs22},numeric);

You have 3 independent variables. You have appearing X(T,Z) and Y(T,R), so 3 independent variables T,Z,R.
That is one two many for pdsolve/numeric.
Did you mean Z when you wrote R?
The syntax would have been
pdsolve({EQ[1],EQ[2]}, BC[1] union BC[2],numeric);

What is Xstar?

1. cosine and sine only appear as cos(g) and sin(g) and g:=Pi/6, so its irrelevant that they appear in input.
2. I don't understand what this is doing in your question: "Failed to load the worksheet /maplenet/convert/Sc.mw . "
3. I looked at your worksheet. Since immediate success is not achieved even after using a midpoint method, as you are told to do by the error message (add method=bvp[midrich] ), and since you are using continuation, I suggest beginning more modestly.
Eq1 has only got the dependent variable f. Thus I suggest studying that first with the boundary conditions
bcsf:=f(0) = 0, D(f)(0) = lambda, D(f)(5) = 1;
#
Continuation in the parameter d will only work if dsolve can solve the problem when d=0, because it works from there and up to d=1 in steps using previously found results as approximate solutions along the way. The problem you really want to solve corresponds to d=1.
So begin by trying d = 0 and one of your lambda values (don't assign to lambda, use eval).
# You can try
res:=dsolve(eval({Eq1,bcsf},{lambda=-2.6,d=0}),numeric,method=bvp[midrich]);
# but you get an error saying

Error, (in dsolve/numeric/bvp) unable to achieve requested accuracy of 0.1e-5 with maximum 128 point mesh (was able to get 0.68e-2), consider increasing `maxmesh` or using larger `abserr`

This, however, may very well lead to a wild goose chase: increasing maxmesh without success. Increasing abserr all the way up to 0.04 works, though:
res:=dsolve(eval({Eq1,bcsf},{lambda=-2.6,d=0}),numeric,method=bvp[midrich],abserr=.04);
##
There is a long way from that to d=1, i.e. using continuation. No success here:
res:=dsolve(eval({Eq1,bcsf},{lambda=-2.6}),numeric,method=bvp[midrich],abserr=.04,continuation=d);
## One may wonder if continuation in the coefficient of eta*f'''(eta) really is a good idea.
##
I tried with (in my view) the absurdly large absolute error abserr=3 (and on the real problem, i.e. d=1):
res1:=dsolve(eval({Eq1,bcsf},{lambda=-2.6,d=1}),numeric,method=bvp[midrich],abserr=3);
#Success, but is it? What I got was
plots:-odeplot(res1,[seq([eta,diff(f(eta),[eta$i])],i=0..2)],legend=[f,seq(Diff(f,[eta$i]),i=1..2)]);


####### Pursuing this further still with abserr=3, no continuation, and lambda=-2.6:
resfull:=dsolve(eval(dsys,{lambda=-2.6,d=1}),numeric,method=bvp[midrich],abserr=3);
plots:-odeplot(resfull,[[eta,f(eta)],[eta,phi(eta)],[eta,theta(eta)]]);


This works for the two other values of lambda as well, but it does bother me that I have to use abserr=3.

First 40 41 42 43 44 45 46 Last Page 42 of 158