Preben Alsholm

13613 Reputation

22 Badges

19 years, 229 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

You can plot and integrate, but I don't see how you can use dsolve with this object.

restart;
M := Array(1 .. 10, 1 .. 2):

for i to 10 do
    M[i, 1] := i;
    M[i, 2] := 3*i;
end do:

M;
Mt := Interpolation:-LinearInterpolation(M);

Mt(3/2);
plot(Mt(t),t=0..10);
int(Mt,0..10,numeric);

 

The problem is that fun is defined outside the procedure. The x and y there are global variables and has nothing to do with the x and y appearing as the parameters in the procedure.
The remedy is simple:
 

restart;
temp_proc := proc(x, y)
local out1, out2, out3, fun:
fun := x^2+y^2;
if x > 0 then out1 := fun; out2 := 2*fun; out3 := k*fun;
elif x <= 0 then out1 := fun; out2 := -2*fun; out3 := -k*fun;
end if:

return(out1, out2, out3);
end proc:

xt := -1: yt := 2:
out1_fin := temp_proc(xt, yt)[1];
out2_fin := temp_proc(xt, yt)[2];
out3_fin := temp_proc(xt, yt)[3];

In the procedure k is treated as global, so out3_fin depends on whether or not k has been assigned a value outside.

maplemint(temp_proc);


 

I wouldn't use DEplot myself. Here is what I would do:

restart;
with(plots):
sys := {diff(x(t), t) = y(t), diff(y(t), t) = -x(t) - 1/2*y(t)}
ics:={x(0) = 1, y(0) = 0};
sol:=dsolve(sys union ics);
solx,soly:=op(subs(sol,[x(t),y(t)]));
plot([solx,soly,t=0..15],labels=[x,y],color=blue);
animate(plot,[[solx,soly,t=0..T],labels=[x,y],color=red,thickness=2],T=0..20);
###################################
## A nonlinear system. We use numeric solution here.
sys2 := {diff(x(t), t) = y(t), diff(y(t), t) = -x(t) - 1/2*y(t)^2}
ics:={x(0) = 1, y(0) = 0};
res:=dsolve(sys2 union ics,numeric);
odeplot(res,[x(t),y(t)],0..200,numpoints=10000,color="DarkGray");
odeplot(res,[x(t),y(t)],0..200,numpoints=10000,color="SkyBlue",frames=50);

MaplePrimes24-04-01_spiral.mw

 

This is a bug.
 

restart;
A := abs(x*Heaviside(x) + x*Heaviside(-x));
plot(A,x=-1..1); #OK
int(A,x=-1..1); # Wrong: should be 1
int(A,x=-1..1,numeric); # 1.000000000
int(A,x=-1..1,method=_RETURNVERBOSE);
int(A,x=-1..1,method=_VERIFYFLOAT);
## Workaround:
B:= x*Heaviside(x) + x*Heaviside(-x);
BP:=convert(B,piecewise);
plot(BP,x=-1..1);
plot(abs(BP),x=-1..1);
int(abs(BP),x=-1..1); # 1

I shall submit a bug report.

In Maple 2024 use adaptive=true, or simply adaptive.
 

res := t^3 - 3*t^2*sqrt(t^2*(12*sqrt(2)*ln(t) + 9*ln(t)^2 + 8)^(1/3))*ln(t)/(12*sqrt(2)*ln(t) + 9*ln(t)^2 + 8)^(2/3) - 2*t^2*sqrt(t^2*(12*sqrt(2)*ln(t) + 9*ln(t)^2 + 8)^(1/3))*sqrt(2)/(12*sqrt(2)*ln(t) + 9*ln(t)^2 + 8)^(2/3);

evalf(eval(res,t=-5));
plot(res,t=-5..1,adaptive,thickness=3);

 

Inspired by Rouben and using his example:
 

restart;
plot([10 + sin(x),0], x=3*Pi..4*Pi);
plot([10 + sin(x),0], x=3*Pi..4*Pi,color=[blue,white],thickness=[3,1]);
plot([-2 + sin(x),0], x=3*Pi..4*Pi);
plot([-2 + sin(x),0], x=3*Pi..4*Pi,color=[blue,white],thickness=[3,1]);
plot([sin(x),0], x=Pi..4*Pi,color=[blue,white],thickness=[3,1]);

 

I'm not sure what you mean by "Kronecker not recognized".
The issue seems to be the following:
 

restart;
A:=Int(
         Sum(p__n*exp(-I*theta)^n, n = 0 .. infinity)*Sum(p__m*exp(theta*I)^m, m = 0 .. infinity), 
       theta = 0 .. 2*Pi);
eval(A,Int=int); # 0 Wrong!
# Counter example:
B:=subs({p__n=exp(-n),p__m=exp(-m)},A); # Don't use eval instead of subs
value(B); #  2*Pi*exp(2)/(exp(2) - 1)

 

A very useful feature of display is overrideoptions:
 

restart;
pt:=plottools:-line([1,2], [3,4],color=red,thickness=2,linestyle=dash); 
plots:-display(pt);
plots:-display(pt,overrideoptions,linestyle=1,thickness=5,color=blue);

restart;
with(IntegrationTools):
J1 := Int(h(x), x=-infinity..a) - Int(h(x), x=-infinity..b);
Ja:=Int(h(x), x=-infinity..a);
Jb:=Int(h(x), x=-infinity..b);
# Case 1: a>b.
Sab:=Split(Ja,b) assuming a>b;
Combine(Sab-Jb);
# Case 2: a<b. 
Sba:=Split(Jb,a) assuming a<b;
Combine(Ja-(Sba)) assuming a<b;
simplify(%);

MaplePrimes24-02-27_IntegrationTools_Combine.mw

NOTE: You don't need the assumption at the end, but it doesn't hurt.
You could add Flip(%); after simplify(%);
Then the results from both cases look the same.

restart;
###
with(LinearAlgebra):
###
v := <x, y>;
w:= <a,b>;
v^%T.v;
v^%T.w;
BilinearForm(v, v) assuming real;
##
## So you could e.g. do this:
M:=module() option package; export BilinearForm;
      BilinearForm:=proc(v::Vector,w::Vector) v^%T.w  end  proc;
end module;  
##             
with(M);
BilinearForm(<x, y>,<a,b>);
BilinearForm(<x, y,z>,<a,b,c>);
##
LinearAlgebra:-BilinearForm(<x, y>,<a,b>);

 

In the first procedure I would do this:
 

restart;
J1 := proc()
  local z:
  z := proc(u) if not u::realcons then return 'procname'(_passed) end if;
  fsolve(sqrt(x)=u, x) end proc:
  evalf(Int(z(u), u=0..1))
end proc:
##############
J1(); # 0.3333333333

In the next I would do likewise:

J2 := proc()
  local z:
   z := proc(q1, q2) if not [q1,q2]::list(realcons) then return 'procname'(_passed) end if;
     exp(
       2*(
         fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = q1, x)
         *
         fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = q2, x)
       )
       -
       fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = sqrt(q1), x)^2
       -
       fsolve(1/2+(1/2)*erf((1/2)*x*sqrt(2)) = sqrt(q2), x)^2
    )
  end proc:
  evalf(Int(z(q1, q2), q1=0..1, q2=0..1))
end proc;

The execution of J2(), however takes an enormous amount of time (infinity?),
but fsolve has to work 4 times inside evalf/int for each value of q1=0..1, q2=0..1.

## Note:  If you comment out evalf(Int(z(q1, q2), q1=0..1, q2=0..1))  and replace it with
z(1/2,1/5) the execution works and returns 0.7300930564.
##
Note 2: I tried using narrower intervals:
evalf(Int(z(q1, q2), q1=0.1..0.9, q2=0.1..0.9))
J2() returned 0.4412675858 in a short time.

You are defining the procedure f to be the same as the procedure F so why not simply start with f:=F  ?

restart;
f := F;
((f@@(-1))@f)(x);
# Let us assume that y is defined this way

y := ((f@@(-1))@g)(x);
# When g is identical to f I would like to get y=x

'y' = eval(y, g=f); # y=x

A concrete example where F=ln:
 

restart;
f := ln;
((f@@(-1))@f)(x);
# Let us assume that y is defined this way

y := ((f@@(-1))@g)(x);
# When g is identical to f I would like to get y=x

'y' = eval(y, g=f); #y=x

 

I corrected a few errors, and here is what I have:
 

restart;
with(plots):
ode1 := diff(f(x), x$3) + 
        f(x)*diff(f(x), x, x) - lambda*diff(f(x), x) - diff(f(x), x)^2 + 2*beta*g(x) - Fr*diff(f(x), x)^2 
        = 0;
ode2 := diff(g(x), x$2) - diff(f(x), x)*g(x) + f(x)*diff(g(x), x) - 2*beta*diff(f(x), x) - lambda*g(x) - Fr*g(x)^2 
        = 0;
# You have 6 boundary conditions, but you can only have 5, unless one of the constants e.g lambda has to be determined.
# Your boundary conditions all have values 0, 
# which means that f(x) = 0, g(x)= 0 for all x on the interval is a (trivial) solution.
# To get a possible nontrivial solution you need to give an approximate solution:
# See ?dsolve,bvp
bcs := f(0) = 0, D(f)(0) = 0, f(5) = 0, D(f)(5) = 0, g(0) = 0, g(5) = 0;
a1 := [lambda = 0.2, beta = 0.2, Fr = 0.1];
####
sys1:=eval({ode1,ode2,bcs[1..5]},a1);# I arbitrarily removed the sixth bcs (g(5)=0).
res1:=dsolve(sys1,numeric);
odeplot(res1,[x,f(x)],thickness=5,view=-0.01..0.01);# As expected the zero solution.

To find a suitable approximate solution is not easy if you have no idea how f and g are expected to look.
You might have an idea?

Here is a proof that y(t) -> infinity as t -> infinity:
MaplePrimes24-02-07_odesys_infinity.mw

Here is the code directly:
 

restart;

sys_ode := diff(x(t), t) = -x(t)^2/(4*Pi*y(t)*(x(t)^2 + y(t)^2)), diff(y(t), t) = y(t)^2/(4*Pi*x(t)*(x(t)^2 + y(t)^2));
ics := x(0) = 1, y(0) = 1:

ode1,ode2:=sys_ode;

ode2/ode1;

ode:=diff(y(x),x)=-(y(x)/x)^3;

sol:=dsolve({ode,y(1)=1});

#From ode1 it follows that x(t) is decreasing from its initial value 1.
# sol shows that x cannot get lower than 1/sqrt(2), which is lower than 1. 
# Since x(t) > 1/sqrt(2), but decreasing we get by replacing x(t) by 1/sqrt(2):
ineq:=rhs(ode2) > subs(x(t)=1/sqrt(2),rhs(ode2));
ode3:=diff(y(t),t)=lhs(ineq);
sol3:=dsolve({ode3,y(0)=1});
limit(rhs(sol3),t=infinity); #infinity
# Thus y(t) from sys_ode will go to infinity too.

####### Graphical illustrations:
plot(rhs(sol),x=1/sqrt(2)..5,color=red,view=0..3); p1:=%: #Saving the plot

Digits:=15;
Sol := dsolve({sys_ode, ics}, {x(t),y(t)}, numeric,abserr=1e-16,relerr=1e-13,events=[[x(t)=5,halt]]);

Sol(100000);
1/sqrt(2.);
plots:-odeplot(Sol,[t,x(t)-1/sqrt(2)],99000..100000,caption=typeset("The difference between x(t) and ",1/sqrt(2)));
#Since x(t) is decreasing from 1 at t=0, you have to go backwards to get near to 1/sqrt(2)
plots:-odeplot(Sol,[x(t),y(t)], t=-43.1..100, color=blue,view=0..3,numpoints=10000);  p2:=%: # Saving the plot

plots:-display(p1,p2);# The last one is on top, so the color is blue
plots:-display(Array([p1,p2]));

 

Here is a simpler version. Let me know, if the result is what you expected:
 

restart;
with(LinearAlgebra):

xA := 1;
yA := 0;
xB := 0;
yB := 0;
xC := 0;
yC := 1;
Mat := Matrix(3, 3, [xA, xB, xC, yA, yB, yC, 1, 1, 1]);

Miv:=Mat^(-1); #Simpler
#Miv := MatrixInverse(Mat);

Miv^%T; #Simpler
#Transpose(Miv);
Miv^%T. <x, y, 1>;

phi := unapply(%,x,y);

for i to 6 do
    B || i := phi(xA || i, yA || i);
end do;

MaplePrimes24-02-05_LinearAlgebra.mw

1 2 3 4 5 6 7 Last Page 2 of 160