Preben Alsholm

13010 Reputation

22 Badges

18 years, 190 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

I removed the question. It clearly appeared to be a repost and was even titled as such. Therefore I deleted the question.

Further questions to a given thread should be given in that thread, not in a new one.

I believe that only the editor of MaplePrimes can restore that question.
Indeed, in your present question you yourself realized that it was a repost of:
https://www.mapleprimes.com/questions/233262-How-Do-I-Solve-Error-in-Dsolvenumericbvpconvertsys

I prefer using dsolve combined with odeplot as shown in the code:

restart;
#k^2*alpha+omega = kao
#K^2*(gamma+alpha)= Kga
f:=(u,z)->z;
g:=(u,z)->kao*u/Kga-beta*u^3/Kga; # You had a wrong sign on the second term
E:=solve({f(u,z)=0,g(u,z)=0},{u,z},explicit);
e1,e2,e3:=op(map(subs,[E],[u,z]));
sys:=diff(u(xi),xi)=f(u(xi),z(xi)), diff(z(xi),xi)=g(u(xi),z(xi));
res:=dsolve({sys,u(0)=u0,z(0)=z0},numeric,parameters=[u0,z0,kao,Kga,beta]);
res(parameters=[0,0.01,.1,.2,1]);
eval([e2,e3],{kao=0.1,Kga=.2,beta=1});
plots:-odeplot(res,[u(xi),z(xi)],0..40,numpoints=10000);
res(parameters=[0,0.1,-.1,-.2,1]);
eval([e2,e3],{kao=-0.1,Kga=-.2,beta=1});
plots:-odeplot(res,[u(xi),z(xi)],-3..3,numpoints=10000); p1:=%:
res(parameters=[0,-0.1,-.1,-.2,1]);
plots:-odeplot(res,[u(xi),z(xi)],-3..3,numpoints=10000,color=blue); p2:=%:
plots:-display(p1,p2);

The first plot:

You may want to have a look at your system etc. again, though.

I corrected several things. Feel free to ask questions about the solution.
Besides technical matters (Maple matters) I found it reasonable to replace de1 by ode1 as seen below.
 

restart;

de1 := (1 - p)*(diff(x(t), t $ 1) - r*x(t)) + p*(diff(x(t), t $ 1) - r*(1 - (x(t) + y(t))/K) + al*x(t)*v(t));

de2 := (1 - p)*(diff(y(t), t $ 1) + m*y(t)+u*y(t)) + p*(diff(y(t), t $ 1) - al*x(t)*v(t) + m*y(t)+u*y(t));

de3 := (1 - p)*(diff(v(t), t $ 1) + eta*v(t)) + p*(diff(v(t), t $ 1) - B*y(t) + eta*v(t));

ibvc := x(0) = 0.005, y(0) = 0.0007, v(0) = 2;


sys1:=eval([de1,de2,de3],p=1);
sys0:=eval~(sys1,[{y=0,v=0},{x=0,v=0},{x=0,y=0}]);
sys_p:=(1-p)*~sys0+~p*~sys1;
ode1,ode2,ode3:=seq(sys_p[j],j=1..3);
ode1; # replaces de1
de1;
ode2; # = de2
de2;
ode3; # = de3
de3;
r:=0.005; al:=0.002; m:=0.0001; K:=0.0138; B:=0.2; eta:=0.006; u:=0.0001;
n := 4;
############First using dsolve directly (numerically)
res:=dsolve({ode1,ode2,ode3,ibvc},numeric,parameters=[p]);
Q:=proc(p,{scene::list:=[t,x(t)],range::range:=0..100}) if not p::realcons then return 'procname(_passed)' end if;
   res(parameters=[p]);
   plots:-odeplot(res,scene,range,_rest)
end proc;
Q(0.5,color=blue);
Q(0.5,scene=[x(t),y(t),v(t)]);
plots:-animate(Q,[p,range=0..50],p=0..1,trace=24);
plots:-animate(Q,[p,range=0..30,scene=[x(t),y(t),v(t)]],p=0..1,trace=24);
############### Now perturbation:
X := unapply(add(g[k](t)*p^k, k = 0 .. n), t);

Y := unapply(add(h[k](t)*p^k, k = 0 .. n), t);

V := unapply(add(i[k](t)*p^k, k = 0 .. n),t);


DE1 := series(eval(ode1, {x = X,y=Y,v=V}), p = 0, n + 1);

DE2 := series(eval(ode2, {x = X,y=Y,v=V}), p = 0, n + 1);

DE3 := series(eval(ode3, {x = X,y=Y,v=V}), p = 0, n + 1);


##### Initial conditions:
M:=eval(([ibvc]), {x(0) = X(0),y(0)=Y(0),v(0)=V(0)});
M0:=eval(M,p=0);
M1,M2,M3,M4:=seq([g[j],h[j],i[j]](0)=~0,j=1..n);  # n=4
ICS:=M0,M1,M2,M3,M4;
###### Now the loop:
g:='g': h:='h': i:='i': # if you do the loop below again you need this. Do it in any case.
for k from 0 to n do

    IBVC1 := select(has,ICS[k+1],g[k]);

    IBVC2 := select(has,ICS[k+1],h[k]);

    IBVC3 := select(has,ICS[k+1],i[k]);

    s1 := dsolve({coeff(DE1, p, k), op(IBVC1)});

    s2 := dsolve({coeff(DE2, p, k), op(IBVC2)});

    s3 := dsolve({coeff(DE3, p, k), op(IBVC3)});

    g[k] := unapply(value(rhs(s1)), t); # value is used because an inert interval may appear

    h[k] := unapply(value(rhs(s2)), t);

    i[k] := unapply(value(rhs(s3)), t);

end do:


'X(t)' = X(t) + O(p^(n + 1));

'Y(t)' = Y(t) + O(p^(n + 1));

'V(t)' = V(t) + O(p^(n + 1));
######## plot tests:
plot(eval(X(t),p=1),t=0..100);
plots:-display(plot(eval(X(t),p=1),t=0..100),Q(1,color=blue,linestyle=3));
Q(1,scene=[t,x(t)-eval(X(t),p=1)],caption="The difference between X(t) and the numerically computed result");
# Do similarly for Y and V

MaplePrimes2021-11-24_perturbation_homotopy.mw

Use unapply instead to define y:
 

restart;
with(inttrans):
y := unapply(invlaplace(exp(-s)/s, s, t),t);
y(t-1);
y(t-1.1);

 

In order to get an answer from your first use of dsolve, where d1,d2,d3,d4, and d5 are not yet assigned, you need assumptions:
 

restart;
eq := diff(Uy(x), x, x) - piecewise(x < d1, 12*F*x/(E*b*h^3), d1 < x and x < d2, 12*((F + F1)*x - F1*d1)/(E*b*h^3), d2 < x and x < d3, 12*((F + F1 + F2)*x - F1*d1 - F2*d2)/(E*b*h^3), d3 < x and x < d4, 12*((F5 + F4 - F)*x + F*L - F5*d5 - F4*d4)/(E*b*h^3), d4 < x and x < d5, 12*((F5 - F)*x + F*L - F5*d5)/(E*b*h^3), 12*F*(L - x)/(E*b*h^3));
sol:=dsolve({eq, Uy(0) = 0, Uy(L) = 0}) assuming d1<d2,d2<d3,d3<d4,d4<d5;
sol:=value(sol) assuming d1<d2,d2<d3,d3<d4,d4<d5;

 

I also use Windows 10. I got the error both times though. No result.
But then I tried replacing integrand with

integrand2:=simplify(integrand);
and I got
ln(x)/x - 1/3*(9*exp(2)*x - 30*x^2 + 18*exp(2) - 236*x - 450)/((-25 + exp(2) - 5*x)*(x + 3)*x)

both times.

 

interface(version);

`Standard Worksheet Interface, Maple 2021.1, Windows 10, May 19 2021 Build ID 1539851`

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1112. The version installed in this computer is 1069 created 2021, September 26, 10:58 hours Pacific Time, found in the directory C:\Users\pkals\maple\toolbox\2021\Physics Updates\lib\`

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

(((-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)

integrand2:=simplify(integrand);

(1/75)*(-75*(x+3)^2*((-(2/5)*x-2)*exp(2)+x^2+10*x+(1/25)*exp(4)+25)*ln(x)+(-120*x^3-1106*x^2-3510*x-4050)*exp(2)+(12*x^2+54*x+81)*exp(4)+225*x^4+3560*x^3+20990*x^2+54000*x+50625)/((x+3)^2*((-(2/5)*x-2)*exp(2)+x^2+10*x+(1/25)*exp(4)+25)*x^2)

print("First time");
int(integrand2,x);

print("second time");
int(integrand2,x);
 

"First time"

ln(x)/x-(1/3)*(9*exp(2)*x-30*x^2+18*exp(2)-236*x-450)/((-25+exp(2)-5*x)*(x+3)*x)

"second time"

ln(x)/x-(1/3)*(9*exp(2)*x-30*x^2+18*exp(2)-236*x-450)/((-25+exp(2)-5*x)*(x+3)*x)

 


 

Download issue_int_nov_11_2021_by_nm.mw

I can add that if I continue after the computation above with integrand instead of integrand2 then I get no errors either time. The answer using integrand has a different form though:
ln(x)/x + 1/x + (15*x^2 + (-4*exp(2) + 356/3)*x - 9*exp(2) + 225)/(x*(x + 3)*(-25 + exp(2) - 5*x))
But the difference between the two answers simplies to 0.

Note: In my original answer I used the two arguments for numboccur in the wrong order. I mistakenly used the same order as used by ListTools:-Occurrences. This was pointed out by Carl Love.
So my answer below has been edited accordingly.
####

Using ListTools:-Occurrences on a list L, a vector V (rtable), and a set all having floating point elements can give quite different results.
I use the example by vv: L:=[1., HFloat(1.), 1., 0., 0.].
The results for the set should of course be expected to differ from the results for the list and the vector.

restart;
#showstat(ListTools:-Occurrences);
L:=[1., HFloat(1.), 1., 0., 0.];
V:=Vector(L);
S:={L[]};
ListTools:-Occurrences~(1.,[L,V,S]);                          #[2, 3, 1]
ListTools:-Occurrences~(HFloat(1.),[L,V,S]);                  #[1, 3, 1]
ListTools:-Occurrences~(1,[L,V,S]);                           #[0, 3, 0]
ListTools:-Occurrences~(1.,[L,V,S],evalb@`=`);                #[3, 2, 2]
ListTools:-Occurrences~(HFloat(1.),[L,V,S],evalb@`=`);        #[3, 2, 2]
ListTools:-Occurrences~(1,[L,V,S],evalb@`=`);                 #[3, 2, 2]
numboccur~([L,V,S],1.);              #[2, 2, 1]
numboccur~([L,V,S],HFloat(1.));      #[1, 1, 1]
numboccur~([L,V,S],1);               #[0, 0, 0]

Here are two ways:
 

restart;
X:=<$0..20>;
Y:=sin~(X);
X1:=<seq(i,i=0..10,0.5)>;
Y1:=cos~(X1);
p:=plot(X,Y); 
p1:=plot(X1,Y1,color=blue); 
plots:-display(p,p1);
## An attractive alternative:
plot([<X|Y>,<X1|Y1>],color=["DarkRed","DarkGreen"],thickness=3,linestyle=[1,3]);

 

You could use subsop:
 

## Notice first the structure:
lprint(expr);
## Physics:-`*`(sigma__y,sigma__x,rho)-Physics:-`*`(sigma__y,rho,sigma__x)
## So
res:=subsop([1,1]=H,[1,2]=1,expr);
lprint(res);
## Result: Physics:-`*`(H,rho)-Physics:-`*`(sigma__y,rho,sigma__x)

 

An alternative to DEplot is odeplot as used here.
I do use DEplot for the direction field though.
MaplePrimes21-10-22_ode_phseplot.mw

One of the images:

plots:-implicitplot(abs(x+2)+abs(y)=3,x=-5..1,y=-3..3);
plots:-implicitplot(abs(y)-abs(x)=2,x=-3..3,y=-5..5);

 

Since this integral appears numerically nontrivial I tried using dsolve/numeric instead of evalf/int.
The result for the real part found below is -0.0513608600887187 (too many digits surely).

MaplePrimes21-10-19_int_by_dsolve.mw
 

restart;

Digits:=15:  #

W := -16*(x - 1/2)*x*(t + 2*w)*(ln((m^2 + w*x*(-1 + x))/m^2)*ln(-1/(w*(t + w)*x^2 - w*(t + w)*x + m^2*t)) + ln(-x*w^2*(-1 + x))*ln((x^2 - x)*w + m^2) + ln(1/m^2)*ln((t + w)*(-1 + x)) + ln(-1/(x*w))*ln(m^2) + dilog(x*w^2*(-1 + x)/(w*(t + w)*x^2 - w*(t + w)*x + m^2*t)) - dilog(x*w*(-1 + x)*(t + w)/(w*(t + w)*x^2 - w*(t + w)*x + m^2*t)))/t^3 + 8*(x - 1/2)*((t + w)*x - t/2)*(-t*(w*x^2 + m^2 - w*x)*(t + w)*ln((x^2 - x)*w + m^2) + w*(x*w*(-1 + x)*(t + w)*ln((t + w)/w) + m^2*t*ln(m^2)))*(2*w*x + t)/(w*(t + w)*(w*(t + w)*x^2 - w*(t + w)*x + m^2*t)*t^3*x) + (16*x^2 - 8*x)/(t*w*(-1 + x)) + (2*t*x + 2*w*x - t - 2*w)*(2*w*x + t - 2*w)*(1 - 2*x)*ln((t*x + (-1 + x)*w)/((-1 + x)*w))/(((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)*t^2) + (-1 + 2*x)*(2*m^2 + w*x - w)*(2*w*x^2 + 2*m^2 - 3*w*x + w)*ln(m^2/(m^2 + w*x*(-1 + x)))/(w^2*(-1 + x)^2*((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)) + (-1 + 2*x)*(2*m^2 + w*x - w)*(2*w*x^2 + 2*m^2 - 3*w*x + w)*ln(m^2/(m^2 + w*x*(-1 + x)))/(w^2*(-1 + x)^2*(w*(t + w)*x^2 - w*(t + w)*x + m^2*t)) + (2*w*x + t)*((t + w)*x - t/2)*(x - 1/2)*ln(w^4/(t + w)^4)/(t^2*(w*(t + w)*x^2 - w*(t + w)*x + m^2*t)) - 8*(x - 1/2)*((-2 + 2*x)*w + t)*((-1 + x)*w + (x - 1/2)*t)*((t*x + (-1 + x)*w)*w^2*(-1 + x)^2*ln((t*x + (-1 + x)*w)/((-1 + x)*w)) + (-(t*x + (-1 + x)*w)*((x^2 - x)*w + m^2)*ln((x^2 - x)*w + m^2) + ln(m^2)*m^2*w*(-1 + x))*t)/((t*x + (-1 + x)*w)*w*(-1 + x)*t^3*((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)) + 16*(x - 1/2)*(ln(1/m^2)*ln((t*x + (-1 + x)*w)*(-1 + x)*w/((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)) + ln((x^2 - x)*w + m^2)*ln(1/((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)) + ln((-1 + x)^2*w^2)*ln((x^2 - x)*w + m^2) - dilog((t*x + (-1 + x)*w)*(-1 + x)*w/((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)) + dilog((-1 + x)^2*w^2/((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)))*(t*x + 2*(-1 + x)*w)/t^3:

Wf:=unapply(W,m,s,t,w): # W as a function of (m,s,t,w), not x

singular(Wf(1,3,-2,-1),0..1);

# With f(x) = Wf(1,3,-2,-1) we split:

IntegrationTools:-Split(Int(f(x),x=0..1),[1/10,1/3,2/3]);

ode:=diff(u(x),x)=Re(Wf(1,3,-2,-1)):

# Using nonsingular initial points x0 = 1/10 and x0 = 2/3

res:=dsolve({ode,u(x0)=0},numeric,parameters=[x0],abserr=1e-15,relerr=5e-10);

# In res Int(f(x),x=x0..x) is computed.

res(parameters=[1/10]);  # setting initial point x0=1/10

v1:=subs(res(1/3),u(x)); # Int(f(x), x = 1/10 .. 1/3)

v2:=subs(res(0),u(x));   # -Int(f(x), x = 0 .. 1/10)

res(parameters=[2/3]);   # setting initial point x0=2/3

v3:=subs(res(1),u(x));   # Int(f(x), x = 2/3 .. 1)

v4:=subs(res(1/3),u(x)); # -Int(f(x), x = 1/3 .. 2/3)

val:=v1-v2+v3-v4;        # -0.0513608600887187

Int(f(x),x=0..1) = val;

 


 

Download MaplePrimes21-10-19_int_by_dsolve.mw

 

It appears that the fact that a is assigned to 1.4 makes a difference to seq. Try using another index:
 

seq([g(i, x[0], y[0])], i = 0 .. 3); # i instead of a

P.S.

The real reason is that the index 'a' also appears in the definition of f:
 

f := (n, m) -> (-n^2 + b*m + a, n);

It is OK that a is assigned as far as seq is concerned, but that 'a' appears in f is the problem.

Since your piecewise expression begins with
 

piecewise(t < t1, kmax, t1 <= t, -kmax, ...)

and since we must either have t < t1 or t1 <= t, the rest of the expression is never used.
P.S. Before thinking about integrating you should try plotting your expression.
Just pick kmax to be e.g. 1 and try different values of t1. Then you will be able to see if you are getting the step function you intended.

A little experimentation:
 

k := piecewise(t < t1, kmax, t1 <= t, -kmax, 2*t1 + sqrt(2)*t1 <= t, kmax, 3*t1 + 2*sqrt(2)*t1 <= t, -kmax, 4*t1 + 2*sqrt(2)*t1 <= t, 0);
plot(eval(k,{kmax=1,t1=2}),t=-10..10);
plots:-animate(plot,[eval(k,{kmax=1}),t=-10..10],t1=-5..5);
## Maybe this is the ordering that you want (?):
k2 := piecewise(2*t1 + sqrt(2)*t1 <= t, kmax, 3*t1 + 2*sqrt(2)*t1 <= t, -kmax, 4*t1 + 2*sqrt(2)*t1 <= t, 0,t < t1, kmax, t1 <= t, -kmax);
plot(eval(k2,{kmax=1,t1=2}),t=-10..10);
plots:-animate(plot,[eval(k2,{kmax=1}),t=-20..20],t1=-5..5,frames=100);

You may also try this:
 

int(eval(k,t1=2), [t = 0 .. t, t = 0 .. t, t = 0 .. t, t = 0 .. t]);
int(eval(k2,t1=2), [t = 0 .. t, t = 0 .. t, t = 0 .. t, t = 0 .. t]);

And finally this:
 

int(k2, [t = 0 .. t, t = 0 .. t, t = 0 .. t, t = 0 .. t]) assuming t1>0,t>0;

 

Would this be OK:

 

Explore(plot(a*x^2+b*x^3,x=-2..2,caption=typeset("a^2+b^2 = ",a^2+b^2)),
parameters = [a = -1.00..1.00, b=-1.00..1.00],placement=right );

 

2 3 4 5 6 7 8 Last Page 4 of 156