Preben Alsholm

MaplePrimes Activity


These are answers submitted by Preben Alsholm

Disregarding your procedure, try this
 

restart;
xr:=1:ao:=sqrt(1+c^2):theta:=arctan(c):a:=ao*exp(I*theta):b:=I*0.5*Delta-a*(k-1)*xr*0.5:no:=1:AA:=5:theta1:=0:Omega:=10:
f:=sqrt(Pi/ao)*exp(-I*0.5*theta)*add(exp(b^2/a)*exp(-a*(k-1)^2*xr^2),k=1..1):alpha:=AA*exp(I*theta1):
f1:= AA^2+((Re(f))^2+(Im(f))^2)*Omega^2+2*Omega*Im(conjugate(alpha)*f):
P1:=plot3d(f1,Delta=-5..5,c=-5..5,axes=boxed,font=[1,1,18]);
##
plottools:-getdata(P1);
A:=%[-1];
mx:=max(A);
tr:=plottools:-transform((x,y,z)->[x,y,z/mx]):
tr(P1);

 

As far as I can see you have exactly the same system as in your previous question:
https://mapleprimes.com/questions/224168-How-To-Solve-Error-Occur-In-DEtools

The difference is only with the intial conditions which are now all zero:
S(0) = 0, V(0) = 0, C(0) = 0, I(0) = 0, R(0) = 0.
As Kitonum points out this obviously gives rise to an error since initially you will be dividing by zero.
Using dsolve/numeric and making I(0) =0.01 gives this (using my code from the link):

ics:={S(0) = 0, V(0) = 0, C(0) = 0, I(0) = 0.01, R(0) = 0};
res:=dsolve(sys union ics, numeric); 
plots:-odeplot(res,[C,I,S](t),3000..3500,numpoints=10000);
plots:-display(Array(1..5,1..1,[seq(plots:-odeplot(res,[t,F(t)],3000..3500,numpoints=10000),F=[C,I,R,S,V])]));

Adjust the time range to your liking. I have 3000..3500.

Here is a way:
 

restart
f:=(3*beta[11]^2-4*beta[11]*sigma[11]+6*beta[12]^2-12*beta[12]*sigma[12]+2*sigma[11]^2+6*sigma[12]^2)*(1/sqrt(6*beta[11]^2-8*beta[11]*sigma[11]+12*beta[12]^2-24*beta[12]*sigma[12]+4*sigma[11]^2+12*sigma[12]^2))*(1/omega^2);
r:=op(indets(f,radical));
subs(r=1/phi,f);

 

Since Kitonum has already dealt with the first system I shall only consider the second:
 

restart;
with(plots):
N := 4: 
de2 := A*(diff(f(eta), eta, eta, eta))+n*(-(diff(f(eta), eta, eta)))^(n-1)*(diff(f(eta), eta, eta, eta))-S*(diff(f(eta), eta))+(2-n)*eta*(diff(f(eta), eta, eta))/(1+n)+2*n*f(eta)*(diff(f(eta), eta, eta))/(1+n)-(diff(f(eta), eta))^2-g(eta)*(diff(f(eta), eta, eta))+(M*M)*(diff(f(eta), eta)) = 0, A*(diff(g(eta), eta, eta, eta))+(-(diff(f(eta), eta, eta)))^(n-1)*(diff(g(eta), eta, eta, eta))-(n-1)*(diff(g(eta), eta, eta))*(diff(f(eta), eta, eta, eta))*(-(diff(f(eta), eta, eta)))^(n-2)-S*(diff(g(eta), eta))+(2-n)*eta*(diff(g(eta), eta, eta))/(1+n)+2*n*f(eta)*(diff(g(eta), eta, eta))/(1+n)-(diff(g(eta), eta))^2+g(eta)*(diff(g(eta), eta, eta))-(M*M)*(diff(g(eta), eta)) = 0, (1+E*j(eta))*(diff(j(eta), eta, eta))+E*(diff(j(eta), eta))^2+2*Pr*n*f(eta)*g(eta)*(diff(j(eta), eta))/(1+n)-Pr*S*(2-n)*eta*(diff(j(eta), eta))/(1+n)+Pr*(Nb*(diff(j(eta), eta))*(diff(h(eta), eta))+Nt*(diff(j(eta), eta))^2)+Pr*lambda*j(eta) = 0, diff(h(eta), eta, eta)+2*Le*Pr*n*f(eta)*g(eta)*(diff(h(eta), eta))/(1+n)-Le*Pr*S*(2-n)*eta*(diff(h(eta), eta))/(1+n)+Nt*(diff(j(eta), eta, eta))/Nb = 0, f(0) = 0, (D(f))(0) = 1, g(0) = 0, (D(g))(0) = alpha, (D(j))(0) = -b*(1-j(0))/(1+E*j(0)), (D(h))(0) = -d*(1-h(0)), (D(f))(N) = 0, (D(g))(N) = 0, j(N) = 0, h(N) = 0: 
####
d2 := subs(alpha = .2, M = .4, A = 1, S = .1, n = .5, Pr = 5, E = 1.5, Nb = .5, Nt = .2, Le = 1, lambda = .2, b = 1.2, d = .5, [de2]): 
####
sol:=dsolve(d2,numeric,approxsoln=[f(eta)=tanh(eta),g(eta)=0.18*tanh(eta),h(eta) = 2/3-eta/6, j(eta) = 1-1/16*eta^2],abserr=1e-3,method=bvp[midrich]);
## Two kinds of plots:
odeplot(sol,[seq([eta,F(eta)],F=[f,g,h,j])]);
display(Array( [seq( odeplot(sol,[eta,F(eta)]),F=[f,g,h,j])]));

The first plot is shown below. MaplePrimes doesn't want an array of plots, so I cannot show the second plot.

What are you going to do with the matrix? Much depends on that!
Wouldn't it be enough to replace wrong values with undefined as in the following simple example:
 

restart;
X:=Vector([$1..10],datatype=float[8]);
Y:=sin~(X);
M:=<X|Y>; #The matrix
plot(M);
M[3,2]:=undefined; # Replace element (3,2) by undefined.
plot(M);
min[defined](M);
max[defined](M);

 

Could you give an example where the need for this comes up?
In fact if f and g are given functions (procedures) then you could just define the composition the obvious way:

p:=x->f(g(x));

If the situation is that f is a function and gx is an expression in x as in the following example then you have a need for unapply:
 

gx:=x^2;
p:=unapply(f(gx),x);

No @ is present here.

@9009134 Is this the correct interpretation of what you got in your worksheet?
 

restart;
CK := .3; Z := 10; L := 1; alpha := .95;
H:=(f,alpha,x,t)->Int((t-s)^(alpha-1)*f(x, s)/GAMMA(alpha), s = 0 .. t); #Changed xi to s
PDE := diff(theta(xi, beta), beta, beta)+L*diff(theta(xi, beta), xi, beta, beta)+diff(theta(xi, beta), beta$3)+(1/2)*diff(theta(xi, beta), beta$4) 
= H(theta,alpha-1,xi,beta)*CK*diff(theta(xi, beta), xi, xi)+H(theta,alpha-1,xi,beta)*(CK*Z+1)*diff(theta(xi, beta), xi, xi, beta)+H(theta,alpha-1,xi,beta)*Z*diff(theta(xi, beta), xi, xi, beta, beta);
Init := {theta(xi, 0) = 0, D[2](theta)(xi, 0) = 0};
Bdry := {theta(0, beta) = 1, theta(10, beta) = 0};
pdsolve(PDE, Init, Bdry, numeric); #Cannot handle the integrals

Notice that I defined H as a function of (f, alpha, x, t) and that in your shorthand version I assumed that f = theta, x = xi, and t = beta.
Anyway, pdsolve cannot solve this type of problems.
##
I'm worried about the denominator in the integrals: (beta-s)^(1.05). This denominator is bound to give problems at s = beta even for smooth theta(xi,s), unless theta(xi,beta)=0 for all xi and beta!
 

There are two problems:

1. The derivatives of lower order are also arguments to D[i](a).
2. diff(u(x,y,t),x,x) = diff(diff(u(x,y,t),x),x). Thus diff(u(x,y,t),x) is an argument of diff.
A solution:
 

restart;
Lambda:=-(1/8)*(D[4](a))(x, y, t, u(x, y, t), diff(u(x, y, t), x), diff(u(x, y, t), y), diff(u(x, y, t), t), diff(u(x, y, t), x, x), diff(u(x, y, t), y, x), diff(u(x, y, t), x, t), diff(u(x, y, t), y, t), diff(u(x, y, t), t, t))*(diff(u(x, y, t), x))-(1/8)*(D[6](a))(x, y, t, u(x, y, t), diff(u(x, y, t), x), diff(u(x, y, t), y), diff(u(x, y, t), t), diff(u(x, y, t), x, x), diff(u(x, y, t), y, x), diff(u(x, y, t), x, t), diff(u(x, y, t), y, t), diff(u(x, y, t), t, t))*(diff(u(x, y, t), y, x))-(1/8)*(D[7](a))(x, y, t, u(x, y, t), diff(u(x, y, t), x), diff(u(x, y, t), y), diff(u(x, y, t), t), diff(u(x, y, t), x, x), diff(u(x, y, t), y, x), diff(u(x, y, t), x, t), diff(u(x, y, t), y, t), diff(u(x, y, t), t, t))*(diff(u(x, y, t), x, t))-(1/8)*(D[5](a))(x, y, t, u(x, y, t), diff(u(x, y, t), x), diff(u(x, y, t), y), diff(u(x, y, t), t), diff(u(x, y, t), x, x), diff(u(x, y, t), y, x), diff(u(x, y, t), x, t), diff(u(x, y, t), y, t), diff(u(x, y, t), t, t))*(diff(u(x, y, t), x, x))-(1/8)*(D[12](a))(x, y, t, u(x, y, t), diff(u(x, y, t), x), diff(u(x, y, t), y), diff(u(x, y, t), t), diff(u(x, y, t), x, x), diff(u(x, y, t), y, x), diff(u(x, y, t), x, t), diff(u(x, y, t), y, t), diff(u(x, y, t), t, t))*(diff(u(x, y, t), x, t, t))-(1/8)*(D[10](a))(x, y, t, u(x, y, t), diff(u(x, y, t), x), diff(u(x, y, t), y), diff(u(x, y, t), t), diff(u(x, y, t), x, x), diff(u(x, y, t), y, x), diff(u(x, y, t), x, t), diff(u(x, y, t), y, t), diff(u(x, y, t), t, t))*(diff(u(x, y, t), x, x, t))-(1/8)*(D[11](a))(x, y, t, u(x, y, t), diff(u(x, y, t), x), diff(u(x, y, t), y), diff(u(x, y, t), t), diff(u(x, y, t), x, x), diff(u(x, y, t), y, x), diff(u(x, y, t), x, t), diff(u(x, y, t), y, t), diff(u(x, y, t), t, t))*(diff(u(x, y, t), y, x, t))-(1/8)*(D[9](a))(x, y, t, u(x, y, t), diff(u(x, y, t), x), diff(u(x, y, t), y), diff(u(x, y, t), t), diff(u(x, y, t), x, x), diff(u(x, y, t), y, x), diff(u(x, y, t), x, t), diff(u(x, y, t), y, t), diff(u(x, y, t), t, t))*(diff(u(x, y, t), y, x, x))+(1/8*(-3*(D[8](a))(x, y, t, u(x, y, t), diff(u(x, y, t), x), diff(u(x, y, t), y), diff(u(x, y, t), t), diff(u(x, y, t), x, x), diff(u(x, y, t), y, x), diff(u(x, y, t), x, t), diff(u(x, y, t), y, t), diff(u(x, y, t), t, t))-4))*(diff(u(x, y, t), x, x, x))-(1/8)*(D[1](a))(x, y, t, u(x, y, t), diff(u(x, y, t), x), diff(u(x, y, t), y), diff(u(x, y, t), t), diff(u(x, y, t), x, x), diff(u(x, y, t), y, x), diff(u(x, y, t), x, t), diff(u(x, y, t), y, t), diff(u(x, y, t), t, t))+(1/4)*(D[8](a))(x, y, t, u(x, y, t), diff(u(x, y, t), x), diff(u(x, y, t), y), diff(u(x, y, t), t), diff(u(x, y, t), x, x), diff(u(x, y, t), y, x), diff(u(x, y, t), x, t), diff(u(x, y, t), y, t), diff(u(x, y, t), t, t))*(diff(u(x, y, t), x, x, x));

L1:=evalindets(Lambda,{seq(specfunc(D[i](a)),i=1..12)},freeze); #First problem handled
L2:=convert(L1,D); # Second problem handled
coeff(L2,D[1](u)(x,y,t));
thaw(%);
coeff(L2,D[1,2](u)(x,y,t));
thaw(%);

An automated version giving the result:
 

S:=indets(L2,function);
for T in S do
  [convert(T,diff),thaw(coeff(L2,T))]
end do;

You can use frontend instead.
 

frontend(coeff,[Lambda,diff(u(x,y,t),x)]);
frontend(coeff,[Lambda,diff(u(x,y,t),x,t)]);
frontend(coeff,[Lambda,diff(u(x,y,t),x,x,x)]);

Here frontend is clearly easiest, but in other situations I have found myself rather frustrated because I couldn't get it to work as I thought it should.

 

Here is yet another solution.
 

restart;
de := diff(z(x),x,x) + z(x) - cos(2*x)/(1+delta*z(x)) = 0; # The given ode
ODE:= diff(z(x),x,x) + z(x) = f(x); # f unassigned
value(dsolve({ODE,z(-Pi/4)=0,z(Pi/4)=0}));
res:=collect(rhs(%),[sin(x),cos(x)],combine); #General formula
# Using:
(1+delta*z(x))^(-1)=convert(1/(1+delta*z(x)),FormalPowerSeries,delta);
F:=x->cos(2*x)*sum((-z(x))^k*delta^k, k = 0 .. infinity);
## The orders 0,1,2 come fast; we must wait a little for order 3:
N:=3: #Don't need to define Z[-1]
for n from 0 to N do
  temp:=eval[recurse](res,{f=F,z=Z[n-1],infinity=n});
  temp:=collect(temp,delta,combine);
  temp:=remove(t->degree(t,delta)>n,temp);
  Z[n]:=unapply(temp,x);
  print(Z[n](x))
end do:

Comparison to the numerically found solution.
 

Resfu:=d->dsolve({eval(de,delta=d),z(-Pi/4)=0,z(Pi/4)=0},numeric);
## The error in using order 3 in delta for delta=0.1:
plots:-odeplot(Resfu(0.1),[x,z(x)-eval(Z[3](x),delta=0.1)],caption=typeset("Error for third order in ",delta));
# Animating in delta from 0.1 and down to 0:
plots:-animate(plots:-odeplot,['Resfu(delta)',[x,z(x)-Z[3](x)]],delta=0.1..0,caption=typeset("Error for third order in ",delta));

The single error plot:

The error for orders 0,1, 2, and 3 when delta = 0.1:
 

plots:-display(Array(
   [seq(plots:-odeplot(Resfu(0.1),[x,z(x)-eval(Z[k](x),delta=0.1)],caption=typeset("Error for order ",k)),k=0..3)]
   ));

 

As you notice simplify isn't always the command that simplifies the most.
In this case, besides factor, radnormal does the job:

radnormal(S);

 

You don't tell us what your answer is as done by hand.
Now do you want to do the same calculation in Maple or do you want to compare your answer to the numerical solution found by Maple?
I shall do the latter and let the solution for delta=0 act as a stand in for your perturbation solution.
 

restart;
ode:=diff(z(x), x, x)-z(x) - cos(2*x)/(1+delta*z(x)) = 0;
res:=d->dsolve({eval(ode,delta=d),z(-Pi/4)=0,z(Pi/4)=0},numeric);
plots:-odeplot(res(0.1),[x,z(x)]);
sol0:=dsolve({eval(ode,delta=0),z(-Pi/4),z(Pi/4)=0}); #result for delta=0
plots:-odeplot(res(0.1),[x,z(x)-rhs(sol0)]);
plots:-animate(plots:-odeplot,['res(delta)',[x,z(x)-rhs(sol0)]],delta=0.1..0);

 

You are probably using Maple 2016 or earlier because in Maple 2017.3 and in the new Maple 2018 I don't get an answer: I get NULL.
Anyway, use value like this:
 

restart;
t20 := r*(diff(theta2(r), r, r))+diff(theta2(r), r)+r*Nb*(-(1/2)*r+(1/4)*h+(1/4)*R0)*((1/16)*Nb*r^3-(1/16)*Nt*r^3-(1/12)*Nb*R0*r^2-(1/12)*Nb*h*r^2+(1/12)*Nt*R0*r^2+(1/12)*Nt*h*r^2+(1/32)*Nb*R0^2*r+(1/16)*Nb*R0*h*r+(1/32)*Nb*h^2*r-(1/32)*Nt*R0^2*r-(1/16)*Nt*R0*h*r-(1/32)*Nt*h^2*r-(1/2)*r-(1/288)*(Nb*R0^4+Nb*R0^3*h-Nb*R0*h^3-Nb*h^4-Nt*R0^4-Nt*R0^3*h+Nt*R0*h^3+Nt*h^4)/((ln(R0)-ln(h))*r)+(1/4)*h+(1/4)*R0)+r*Nb*(-(1/4)*R0+(1/2)*r-(1/4)*h+(1/4)*Nt*R0/Nb-(1/2)*Nt*r/Nb+(1/4)*Nt*h/Nb)*(-(1/4)*R0+(1/2)*r-(1/4)*h)+2*r*Nt*(-(1/4)*R0+(1/2)*r-(1/4)*h)*((1/16)*Nb*r^3-(1/16)*Nt*r^3-(1/12)*Nb*R0*r^2-(1/12)*Nb*h*r^2+(1/12)*Nt*R0*r^2+(1/12)*Nt*h*r^2+(1/32)*Nb*R0^2*r+(1/16)*Nb*R0*h*r+(1/32)*Nb*h^2*r-(1/32)*Nt*R0^2*r-(1/16)*Nt*R0*h*r-(1/32)*Nt*h^2*r-(1/2)*r-(1/288)*(Nb*R0^4+Nb*R0^3*h-Nb*R0*h^3-Nb*h^4-Nt*R0^4-Nt*R0^3*h+Nt*R0*h^3+Nt*h^4)/((ln(R0)-ln(h))*r)+(1/4)*h+(1/4)*R0);
t22 := theta2(h) = 0, theta2(R0) = 0;
res:=dsolve({t20,t22});
resv:=value(res) assuming r>h,h>0,R0>h;
odetest(resv,[t20,t22]); #Check OK

 

In your procedure replace nops by numelems.
 

restart;
#Example:
v:= <1, 2, 3, 2, 4>;
nops(v);
numelems(v);
op(v);
for i from 1 to nops(v) do op(i,v) end do;

 

I shall assume that you only want the answer for sin(x) >=0.
Thus I shall only consider 0<= T <= Pi.
 

restart;
res:=int(sqrt(sin(x)),x=0..T) assuming T>0, T<Pi;
diff(res,T);
simplify(%) assuming T>0,T<Pi; #WRONG on Pi/2..Pi
## Thus we split:
res1:=int(sqrt(sin(x)),x=0..T) assuming T>0, T<Pi/2;
simplify(diff(res1,T)) assuming T>0,T<Pi/2; #OK
res2:=int(sqrt(sin(x)),x=Pi/2..T) assuming T>Pi/2,T<Pi;
simplify(diff(res2,T)) assuming T>Pi/2,T<Pi; #OK
resp:=piecewise(T<Pi/2,res1,eval(res1,T=Pi/2)+res2); # Result on 0..Pi
plot(resp,T=0..Pi);

The title of your question was: "More than one solution for a special case."
But fsolve will only give you one result unless the input is a polynomial in one variable.
###
Well, since dsolve gives you a result, try it again with a higher setting of Digits.
Just for a short experiment I set Digits to 30 (you can try much higher):
Removing the imaginary solutions I got these 4 solutions (when B=10):
[{C1 = -8.21367014680801194633270213292*10^(-12), Nu = -4.38002192864507929874554229402*10^6}, {C1 = 20.4996314736857070747869876646, Nu = -65.8211157946078643801270593845}, {C1 = -.898697538127575625586802649353, Nu = -53.0523912091922304391249719870}, {C1 = -8.21426028148815833335058821461*10^(-12), Nu = 4.37969614250226629250101130330*10^6}].
###########
Trying your loop using dsolve instead of fsolve and Digits=50 and also removing imaginary results.
You might as well be using solve. The use of dsolve here is confusing (was to me!).

Bvals:=[ -10.0, -1.0, -0.5, 0, 0.5, 1.0, 10.0
       ]:
for j from 1 to nops(Bvals) do
   temp:= [evalf(dsolve(eval({E1, E2},B=Bvals[j]), indets(eval({E1, E2},B=Bvals[j]))))];
   res[j]:=remove(has,temp,I)
end do:
for j from 1 to nops(Bvals) do nops(res[j]) end do;

I'm getting 2 real results for the first 4 B values and 4 for the last 3 (edited).
I would suggest using solve instead of dsolve.

First 11 12 13 14 15 16 17 Last Page 13 of 149