Preben Alsholm

13663 Reputation

22 Badges

19 years, 343 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@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

@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;

diff(diff(diff(f(eta), eta), eta), eta)+f(eta)*(diff(diff(f(eta), eta), eta))-(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;

diff(diff(theta(eta), eta), eta)+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;

diff(diff(phi(eta), eta), eta)+Le*(diff(phi(eta), eta))+Nt*(diff(diff(theta(eta), eta), eta))/Nb

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

f(0) = 0, (D(f))(0) = 1, (D(f))(10) = 0

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

theta(10) = 0, (D(theta))(0) = -Bi*(1-theta(0))

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

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:

[.1, .2, .3, .4, .5]

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

M := Matrix(5, 6, {(1, 1) = .1, (1, 2) = HFloat(0.09079204195997714), (1, 3) = HFloat(0.07999458907826575), (1, 4) = HFloat(0.06038735511214605), (1, 5) = HFloat(0.035939234854121746), (1, 6) = HFloat(0.017220024220039862), (2, 1) = .2, (2, 2) = HFloat(0.09048653015958226), (2, 3) = HFloat(0.07866376250695364), (2, 4) = HFloat(0.056451951494423856), (2, 5) = HFloat(0.0306001744671172), (2, 6) = HFloat(0.013603325725012235), (3, 1) = .3, (3, 2) = HFloat(0.09015411370183106), (3, 3) = HFloat(0.07711548838707), (3, 4) = HFloat(0.051978605826939764), (3, 5) = HFloat(0.02573745440282669), (3, 6) = HFloat(0.010872451216143206), (4, 1) = .4, (4, 2) = HFloat(0.08979077608034479), (4, 3) = HFloat(0.07529813593135734), (4, 4) = HFloat(0.04711421671313636), (4, 5) = HFloat(0.021624435358676336), (4, 6) = HFloat(0.00885989600363249), (5, 1) = .5, (5, 2) = HFloat(0.08939164233960008), (5, 3) = HFloat(0.07314965034262784), (5, 4) = HFloat(0.04215277835976381), (5, 5) = HFloat(0.018312484329286654), (5, 6) = HFloat(0.0073775356466370765)})

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

top := Vector[row](6, {(1) = Nt, (2) = Nb = .1, (3) = Nb = .2, (4) = Nb = .3, (5) = Nb = .4, (6) = Nb = .5})

Vals:=<top,M>;

Vals := Matrix(6, 6, {(1, 1) = Nt, (1, 2) = Nb = .1, (1, 3) = Nb = .2, (1, 4) = Nb = .3, (1, 5) = Nb = .4, (1, 6) = Nb = .5, (2, 1) = .1, (2, 2) = HFloat(0.09079204195997714), (2, 3) = HFloat(0.07999458907826575), (2, 4) = HFloat(0.06038735511214605), (2, 5) = HFloat(0.035939234854121746), (2, 6) = HFloat(0.017220024220039862), (3, 1) = .2, (3, 2) = HFloat(0.09048653015958226), (3, 3) = HFloat(0.07866376250695364), (3, 4) = HFloat(0.056451951494423856), (3, 5) = HFloat(0.0306001744671172), (3, 6) = HFloat(0.013603325725012235), (4, 1) = .3, (4, 2) = HFloat(0.09015411370183106), (4, 3) = HFloat(0.07711548838707), (4, 4) = HFloat(0.051978605826939764), (4, 5) = HFloat(0.02573745440282669), (4, 6) = HFloat(0.010872451216143206), (5, 1) = .4, (5, 2) = HFloat(0.08979077608034479), (5, 3) = HFloat(0.07529813593135734), (5, 4) = HFloat(0.04711421671313636), (5, 5) = HFloat(0.021624435358676336), (5, 6) = HFloat(0.00885989600363249), (6, 1) = .5, (6, 2) = HFloat(0.08939164233960008), (6, 3) = HFloat(0.07314965034262784), (6, 4) = HFloat(0.04215277835976381), (6, 5) = HFloat(0.018312484329286654), (6, 6) = HFloat(0.0073775356466370765)})

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

[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)
             ]));

 

 

 

 

 

 

 

``

Download MaplePrimes21-12-24_odebvpII_executed.mw

@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))])]);

@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)))
             ]));

I deleted your recent repost of this question. It appears that you already received an answer.
If you have further questions then continue this thread. Don't start another one.

We may use that the system has a conserved quantity.
Also it's worth noticing that the system can de solved exactly. We shall not use that solution though.

## sys is defined above in my answer.
dsolve({sys});
ode:=eval(diff(sys[1],xi),sys[2]);
EQ1:=map(int,ode*diff(u(xi),xi),xi)+(0=C);
isolate(EQ1,C);
EQ2:=subs(sys[1],u(xi)=u,z(xi)=z,%);
params1:={kao=0.1,Kga=.2,beta=1};
plots:-display(seq(plots:-implicitplot(eval(EQ2,params1 union {C=k/200}),u=-1..1,z=-1..1),k=-2..2));
params2:={kao=-0.1,Kga=-.2,beta=1};
plots:-display(seq(plots:-implicitplot(eval(EQ2,params2 union {C=k/10}),u=-1..1,z=-1..1),k=-2..2));
rhs(EQ2);
plots:-contourplot(eval(rhs(EQ2),params1),u=-1..1,z=-1..1,contours=[seq(k/10,k=-2..2)],grid=[50,50]);
plots:-contourplot(eval(rhs(EQ2),params2),u=-1..1,z=-1..1,contours=[seq(k/10,k=-2..2)],grid=[50,50]);

The plots from contourplot are not as good as the ones from implicitplot for C = 0 and (u,z) near (0,0).

MaplePrimes (right now anyway) seems to have a serious problem bringing the (first) implicitplot when I just copy and paste.
So I saved the plot as a png file:


Here is the second implicitplot:

@SuganyaG I knew that my worksheet in my answer above had lots of details that needed explanation.

Here is an attempt:
MaplePrimes2021-11-25_perturbation_2.mw

Comment: I also tried setting n:=5 instead of n:=4.
This works too, but plot gets a problem for some reason, probably because of the huge size of the input.
A simple way out is to set the environmental variable UseHardwareFloats to false, like this
UseHardwareFloats:=false;
Its default value is 'deduced'.

@Carl Love You are quite right. Thank you.

@Carl Love Thank you Carl for pointing out my mistake. I will revise my answer accordingly. This also means that my statement about vectors will be removed.

@nm You are right. Something strange is going on.
If in a new worksheet I start with my integrand2 and then do int(integrand2,x) I get an error.
If I then execute simplify(integrand2) and then int(integrand2,x) I get the result.
And if again in a fresh worksheet I start with your original integrand, then do simplify(integrand) then finally
int(integrand,x) then I get a result. Notice that I'm not using the simplified version:

restart;
integrand:=(((-3*x^2-18*x-27)*exp(2)^2+(30*x^3+330*x^2+1170*x+1350)*exp(2)-75*x^4-1200*x^3-7050*x^2-18000*x-16875)*ln(x)+(12*x^2+54*x+81)*exp(2)^2+(-120*x^3-1106*x^2-3510*x-4050)*exp(2)+225*x^4+3560*x^3+20990*x^2+54000*x+50625)/((3*x^4+18*x^3+27*x^2)*exp(2)^2+(-30*x^5-330*x^4-1170*x^3-1350*x^2)*exp(2)+75*x^6+1200*x^5+7050*x^4+18000*x^3+16875*x^2);
simplify(integrand);
int(integrand,x);

I get the errors also from your Calendar version though.
Using restart doesn't remove everything (see ?restart).
So in testing this we must start with a fresh worksheet (a new engine).

 

@ You are wrong. The general term tends to 0:
 

limit(1/k^(2-cos(1/k)),k=infinity); # 0
## Obvious also since the term is positive and satisfies:
1/k^(2-cos(1/k))<= 1/k; 

 

@Kitonum If the result isn't obvious it may be a good idea to see how Maple finds that:
 

restart;
infolevel[int]:=5;
int(1/k^(2-cos(1/k)), k=1..infinity);
# So try this
asympt(1/k^(2-cos(1/k)),k);

So the result infinity appears to be OK.

@dearcia You must have an older Maple version.
Replace DEtools:-DEplot by DEtools[DEplot] and you should be OK.
If, however, your version is older than Maple 13 then elementwise application (using the tilde ~) is not available. You can work around that rather easily by using map.
###############################
Here is a version that I made for Maple 12 on the basis of my Maple 2021 version.
That made me realize how convenient the elementwise operation is. It came with Maple 13.
There are several changes as you will notice.
## This version also works in Maple 2021.
MaplePrimes21-10-22_ode_phaseplot_VersionM12.mw
 

restart;

F0 := 0; zeta := .25; w := 1; Omega := 1; m := 1;

ode1 := diff(x(t), t) = y(t); ode2 := diff(y(t), t) = F0*cos(Omega*t)/m-2*zeta*w*y(t)-w^2*x(t);

res:=dsolve({ode1, ode2,x(0)=x0,y(0)=y0},numeric,parameters=[x0,y0]);

P:=proc(xy0::list,tR::range:=0..10) if not xy0::list(realcons) then return 'procname(_passed)' end if;
   res(parameters=xy0);
   plots:-odeplot(res,[x(t),y(t)],tR,_rest)
end proc;

#Test of P

P([0,50],linestyle=dash,color=blue,thickness=3);

P([0,50],-5..15,linestyle=dash,color=blue,thickness=3);

L:=[[x(0) = 0, y(0) = 50], [x(0) = 9, y(0) = 25], [x(0) = 85, y(0) = 20], [x(0) = .25, y(0) = .5], [x(0) = 7, y(0) = 5]];

L1:=map(x->map(rhs,x),L);

opts:=zip(`=`,[color$5],[red, green, black, navy, maroon]),[(thickness=3)$5],zip(`=`,[linestyle$5],[dash, dashdot, dot, longdash, spacedash]);

nopts:=nops([opts]);

p0:=DEtools[DEplot]([ode1, ode2], [x(t), y(t)], t = 0 .. 1, x = -100 .. 100, y = -100 .. 100,
 [[x(0) = 0, y(0) = 0]],arrows = medium,color=green):

plots[display](seq(P(L1[i],seq(opts[j,i],j=1..nopts)),i=1..5),p0);

plots[display](seq(P(L1[i],-0.9..10,seq(opts[j,i],j=1..nopts)),i=1..5),p0);

 


 

Download MaplePrimes21-10-22_ode_phaseplot_VersionM12.mw

 

@tomleslie OK, I misread the result from A3 as having a minus Float(infinity).

But my real mistake was to use a wrong ordering of the arguments of f3.
I should have considered
 

W1:=eval(W,[m,s,t,w]=~[1, 3, -2, -1]);
plot(Re(W1),x=0..1);

instead of

W1:=eval(W,[s,t,w,m]=~[1, 3, -2, -1]);

which gives no problems.
The plot of the correct W1 shows a clear indication of a singularity at x = 1/3
 

singular(W1,0..1); # {x = 0}, {x = 1/3}, {x = 1}
evalf(Int(Re(W1),x=0..1/3)); # -0.002885966034
evalf(Int(Re(W1),x=1/3..1)); # Problems

 

@tomleslie Obviously the real value is wrong since it must be larger than -2 if the following plots can be believed:
 

int(-1/(1-x)^(1/2),x=0..1);
plot([Re(eval(W,{s=1,t=3,w=-2,m=-1})),-1/(1-x)^(1/2)],x=0..1);

 

First 17 18 19 20 21 22 23 Last Page 19 of 229