13613 Reputation

19 years, 226 days

D and i...

The use of D is unfortunate since it means the differentiation operator in Maple:

`D(sin); # Answer cos`

So use another name, e.g. D1.
Secondly, by default the imaginary unit is  'I', not 'i'.

An observation...

@acer Just an observation: So the following two versions of dividing are not quite the same:

```1/4/(x-2);
1/(4*(x-2));
```

Correction...

@tomleslie In order to be fair you must divide by 2000 000 in ans2:

```ans1:=(654321.987*123456.789)/2;
ans2:=(654321987*123456789)/2000000;```

Floats and syntax...

@nm Normally you cannot put a number first followed by a letter; that would produre a syntax error. Example 1w would produce an error.
There are exceptions:  1d and  1e are examples.

@Jaime_mc2 This works with Digits:=50 in Maple 2021.2

```evalf(Int(abs(u(x) - P__2), x = -1 .. 1,method= _Gquad));
```

The result is 0.00062767200590347019035308761635491656933673486658951.
Inspired by mmcdara I added simplify like this:

`r := evalf(allvalues(RootOf(simplify(ChebyshevT(5, x)), x)));`

and this time the result was
0.00062767200590347019035308761635491656933673486658622

not a huge difference there.

But the latter version, but not the former, worked in Maple 12 too and with that same result.

The two versions don't work the same way...

@nm Have a look at this simple example:

```restart;
f:=x->a*x;
simplify(sqrt(f(a)),assume=[a<0]);
simplify(sqrt(f(a))) assuming a<0;
```

In the first case sqrt(f(a)) is evaluated to sqrt(a^2) before the assumption kicks in, thus the result -a.

In the second case the expression sqrt(f(a)) is not evaluated until the assumption a<0 is made. But the "version" of the 'a'  that is seen by `assuming` is therefore only the argument to f, not the 'a' appearing in the in the definition of f.
My point here has just been to say that the two versions don't work quite the same.
###
You can see it here from the printout from debug(`assuming`:

```restart;
f:=x->a*x;
debug(`assuming`);
simplify(sqrt(f(a)),assume=[a<0]);
simplify(sqrt(f(a))) assuming a<0;
```

Replace Int with F...

Obviously the syntax for Int is wrong.
Replace Int with F then both simplifications work.

Maple 18...

But in Maple 18 (which you are using) I get the plot you show.

But as Carl Love points out you should use a finer grid. Besides that it helps to set gridrefine to 2 as in this code:

```plots:-display(seq(plots:-implicitplot(eval(EQ2,params1 union {C=k/200}),u=-1..1,z=-1..1,grid=[50,50],gridrefine=2),k=-2..2));
```

If you look at the help for implicitplot you will find other options. Actually it looks to me that this version gives a better result than the one above:

```plots:-display(seq(plots:-implicitplot(eval(EQ2,params1 union {C=k/200}),u=-1..1,z=-1..1,grid=[50,50],gridrefine=1,crossingrefine=2),k=-2..2));
```

Merry Christmas...

@jennierubyjane Merry Christmas to you too!

Not quite right...

@jennierubyjane You can say considerably more.
A few hints at what is going on:
Notice in my code the use of values for initmesh and maxmesh.
So from the interval, in your case 0..10, is chosen an initial mesh (default consisting of 8 points). I used 512 points.
Notice also that an approximate solution can be given as an option (I didn't use that here). If not given the procedure finds an approximate solution itself based only on the boundary conditions (!).
Here is an example of getting some information about the solving process:

```restart;
#infolevel[`dsolve/numeric/BVP`]:=2:  # Also try 3 or 4:
infolevel[`numeric/bvp`]:=5: # Another possibility
ode:= diff(y(x),x,x)+sin(y(x)) +y(x) = 1;
bcs:=y(0)=1,D(y)(5)=0;
res:=dsolve({ode,bcs},numeric);
plots:-odeplot(res,[x,y(x)]);
```

Deliberately I used a much simpler example than yours, but the process is the same.
Notice that `dsolve/numeric/bvp` is different from `dsolve/numeric/BVP`.
The former calls the latter, i.e. `dsolve/numeric/BVP` is used by `dsolve/numeric/bvp`.
You can see that here in line 440.
showstat(`dsolve/numeric/bvp`);

dsolve/numeric/bvp...

@jennierubyjane When dsolve/numeric is given a boundary value problem it is using the procedure `dsolve/numeric/bvp`

For that procedure see the help page ?dsolve,numeric,bvp.

RKF45 is made for initial value problems.
If one insists on somehow making use of that method (or other intial value procedures) shooting has to be done.
That means starting at one end (e.g. eta = 0) with all the required initial values, which implies "guessing" some.

That can be done. I have no reason to believe that that is going to give a better result.
The "guessing" can be implemented by using the parameters option to dsolve where the parameters will be the unknown initial values. Then fsolve will be involved after that.

Revised worksheet...

@jennierubyjane I revised the worksheet so phi and -diff(theta(eta),eta) are shown and for all values also of Nb.

Why the numbers in your table don't agree with the Maple results I cannot tell.
Are the original equations OK?
Was the right endpoint i.e. eta = 10 also used by the author of the table?
I understand that here it stands for infinity. Sometimes I have seen lower values than 10 used, e.g. 8.
The main thing, however, is that we compare with the same values of input all over.
Does the author reveal what kind of method was used to find the values?

MaplePrimes21-12-24_odebvpIII.mw

The full table and colors...

@jennierubyjane I have made several changes as you will notice.

I usually don't upload an executed version of a worksheet, and below you will see that the quality isn't nearly as good as it is on my computer.
When you click on the display/insequence version a panel on top of the screen comes up (or should come up). You can choose the speed or move one frame at a time.
I notice that the display of Vals is screwed up. Have a look at the worksheet itself.

MaplePrimes21-12-24_odebvpII_executed.mw

 > restart;
 > with(plots):
 >
 > DE1 := diff(f(eta), eta\$3)+f(eta)*(diff(f(eta), eta\$2))-(diff(f(eta), eta))^2 = 0;

 > DE2 := diff(theta(eta), eta\$2)+Pr*f(eta)*(diff(theta(eta), eta))+Pr*Nb*(diff(theta(eta), eta))*(diff(phi(eta), eta))+Pr*Nt*diff(theta(eta), eta)^2 = 0;

 > DE3 := diff(phi(eta), eta\$2)+Le*(diff(phi(eta), eta))+Nt*(diff(theta(eta), eta\$2))/Nb;

 > BC1 := f(0) = 0, (D(f))(0) = 1, (D(f))(10) = 0;

 > BC2 := theta(10) = 0, (D(theta))(0) = -Bi*(1-theta(0));

 > BC3 := phi(0) = 1, phi(10) = 0;

 > Digits:=15:
 > ##
 > L := [0.1,0.2,0.3,0.4,0.5]; nL:=numelems(L): for i from 1 to nL do    for k from 1 to nL do       res[i,k]:=dsolve(                   eval({BC1, BC2, BC3, DE1, DE2, DE3}, [Nt = L[k], Nb=L[i], Pr=10, Le=10, Bi=0.1]),                   numeric,output=listprocedure, abserr=1e-10,maxmesh=1024,initmesh=512                 );       PH[i,k],mTH1[i,k] := op(subs(res[i,k], [phi(eta),-diff(theta(eta), eta)]));        end do; end do:

 > M:=Matrix([seq([L[k],seq(mTH1[i,k](0),i=1..nL)],k=1..nL)]);

 > top:=Vector[row]([Nt,seq(Nb=L[i],i=1..nL)]);

 > Vals:=;

 > colorList:=[red,blue,"Green",magenta,"DarkGreen"];

 > display(Array([                display(seq(odeplot(res[1,k],[eta,phi(eta)],0..2.5,color=colorList[k]),k=1..nL)),                          display(seq(odeplot(res[1,k],[eta,theta(eta)],0..2.5,color=colorList[k]),k=1..nL))              ]));

 > display(Array([                display(seq(odeplot(res[1,k],[eta,phi(eta)],0..2.5,color=colorList[k]),k=1..nL),insequence),                          display(seq(odeplot(res[1,k],[eta,theta(eta)],0..2.5,color=colorList[k]),k=1..nL),insequence)              ]));

 >
 >
 >

Range and values...

@jennierubyjane To get another range than the default taken from the boundary conditions, just add a range:

```display(Array([
display(seq(odeplot(res[k],[eta,phi(eta)],0..2.5),k=1..numelems(L))),
display(seq(odeplot(res[k],[eta,theta(eta)],0..2.5),k=1..numelems(L)))
]));
```

To get values basically you can just do as here:

`res[1](0.5); # eta=0.5`

If you only want phi and theta and have a particular vector of eta values in mind then do:

```for k from 1 to numelems(L) do
PH[k],TH[k]:=op(subs(res[k],[phi(eta),theta(eta)]));
end do:
interface(rtablesize=11); # So you actually see all the values (they are there anyway)
ET:=Vector([seq(0.25*i,i=0..10)]);
Matrix([ET,PH[3]~(ET),TH[3]~(ET)]); # Example: k = 3, i.e. Pr = L[3] = 5, Columns: eta,phi,theta
```

If you want one particular value of eta, here eta = 0.5, and phi and eta for all the values of Pr do:

```## Values for all Pr and for eta = 0.5: Pr, phi, theta
Matrix([Vector(L),Vector([seq(PH[k](0.5),k=1..numelems(L))]),Vector([seq(TH[k](0.5),k=1..numelems(L))])]);
```

Values...

@jennierubyjane Which values do you want to print out?
You said:  " ... dont know how to print out the values tsk tsk"
What does "tsk tsk" refer to?

From your do loop in the code I see that the procedures for phi(eta) and -diff(phi(eta),eta) are saved.
Are those numerical procedures the ones you call values?

While waiting for your answer to my question I shall allow myself to continue mmcdara's worksheet
ODEprobz_mmcdara.mw here:

```L := [0.7e-1, .2, .7, 2, 7, 20, 70];
for k from 1 to numelems(L) do
res[k]:=dsolve(
eval({BC1, BC2, BC3, DE1, DE2, DE3}, [Pr = L[k], Nb=1, Nt=1, Le=1, Bi=1]),
numeric,output=listprocedure
);

end do:
display(Array([
display(seq(odeplot(res[k],[eta,phi(eta)]),k=1..numelems(L))),
display(seq(odeplot(res[k],[eta,-diff(phi(eta),eta)]),k=1..numelems(L)))
]));

```

 First 16 17 18 19 20 21 22 Last Page 18 of 229
﻿