Kitonum

15529 Reputation

24 Badges

12 years, 140 days

MaplePrimes Activity


These are answers submitted by Kitonum

Index:=proc(v::{Vector,Array})
local i0,i1,i;
i0:=`if`(v::Vector,1,lhs(op(v)[1]));
i1:=`if`(v::Vector,op(v)[1],rhs(op(v)[1]));
for i from i0 to i1 do
if v[i]>=0 then return i fi;
od;
end proc: 

Examples of use:
V:=Vector([-1,-1,-2,0,3,4]):
Index(V);
V1:=Array(0..3, [-1,-2,4,-5]):
Index(V1);
V2:=Vector([-1,-1,-2]):
Index(V2);

                                 4
                                 2


Edit.
                                       

f:=x->-(x^3-6*x^2+11*x-6)*(x^2-4*x+5)/(x^4-8*x^3+24*x^2-32*x+17);
plot([f(x),2-x], x=-2..6, linestyle=[1,3],color=[red,blue], thickness=[2,0], scaling=constrained);

            


Addition. There is a rational function with polynomials of lower degrees satisfying all conditions:

f:=x->(-x^3+6*x^2-11*x+6)/(x^2-4*x+5);
convert(f(x),parfrac);
plot([f(x),2-x], x=-4..8, -5..5,  linestyle=[1,3],color=[red,blue], thickness=[2,0], scaling=constrained);

 

Here is a procedure for this:

Grouping:=proc(L::listlist)
local L1;
uses ListTools;
L1:=[Categorize((x,y)->x[2]=y[2], L)];
map(p->[p[1,2],map(l->l[1],p)], L1);
end proc:


Example of use:

MyNestedList := [   [g1, [1,2,3]]     ,   [g2, [0,1,3]]   ,  [g3, [1,2,3]]  ,  [g4, [0,1,3]]]:
Grouping(MyNestedList);
                 
 [[[1, 2, 3], [g1, g3]], [[0, 1, 3], [g2, g4]]]

I replaced  D  with  Daw  because  D  is the differential operator in Maple. We have a transcendental equation  m = 0.9058763835e-3*sqrt(6)/(x^2*m^2*exp(-26.00000000*x*m))-1/2-0.2884615385e-1/(x*m)  that we can solve only numerically. On the plot we see that for each  x   there are 2 values of  . In the code we find  m  with the smallest modulus:
 

restart;
Daw:=x->exp(-x)*int(_t^2,_t=0..sqrt(x));
N:=q->exp(q)/sqrt(q)*sqrt(3/2)*Daw(3*q/2);

proc (x) options operator, arrow; exp(-x)*(int(_t^2, _t = 0 .. sqrt(x))) end proc

 

proc (q) options operator, arrow; exp(q)*sqrt(3/2)*Daw((3/2)*q)/sqrt(q) end proc

(1)

t:=0.6: n:=6: k:=1.3:
q:=(24/(n-3))*(x*k*m)/t;
Eq:=m= 1/2*(exp(q)/(q*N(q))-1-1/q);
plots:-implicitplot(Eq, x=0..10,m=-1..0, gridrefine=3);

17.33333333*x*m

 

m = 0.9058763835e-3*6^(1/2)/(x^2*m^2*exp(-26.00000000*x*m))-1/2-0.2884615385e-1/(x*m)

 

 

f:=unapply('fsolve'(m= 1/2*(exp(q)/(q*N(q))-1-1/q),m=-0.25..0),x);

proc (x) options operator, arrow; fsolve(m = 0.9058763835e-3*6^(1/2)/(x^2*m^2*exp(-26.00000000*x*m))-1/2-0.2884615385e-1/(x*m), m = -.25 .. 0) end proc

(2)

# Examples of use
f(0.5), f(1), f(5);

-.2066314659, -0.7733721921e-1, -0.1398887471e-1

(3)

 


 

Download Equation.mw

restart;
F:=()->RandomTools:-Generate(choose([1$2, 2$5, 3, 4, 5])):
# Example of use
seq(F(), i=1..100);
ListTools:-Collect([%]);

4, 1, 4, 2, 1, 2, 1, 1, 4, 2, 2, 2, 5, 1, 2, 5, 2, 3, 2, 2, 2, 2, 

  1, 1, 2, 4, 2, 5, 1, 1, 2, 3, 5, 2, 1, 2, 2, 4, 1, 1, 1, 1, 4, 

  2, 2, 5, 2, 2, 1, 3, 2, 1, 2, 2, 2, 5, 2, 4, 2, 5, 2, 1, 2, 2, 

  4, 1, 5, 4, 2, 1, 3, 2, 3, 3, 5, 2, 2, 1, 2, 1, 2, 1, 2, 2, 3, 

  2, 2, 1, 5, 3, 5, 5, 5, 2, 2, 1, 2, 4, 2, 1
          [[1, 25], [2, 44], [3, 8], [4, 10], [5, 13]]
 

See this a toy example:

A:=a+b*x=0:
coeff(A, x);
coeff(lhs(A), x);

                       Error, unable to compute coeff
                                               b

 

Use  fsolve  with  initialpoint  option :

fsolve([Asubs/Bsubs = 1, 2*Esubs = h*`J__ω`[49]], {a = 2*10^(-20) .. 4*10^(-20), k = 4*10^8 .. 8*10^8});

                            {a = 3.31129407052*10^(-20), k = 6.47179093038*10^8}


I got this this initialpoint  {a = 2*10^(-20) .. 4*10^(-20), k = 4*10^8 .. 8*10^8}  from the implicitplot. If you want to automate this process for a large number of applications, then try using the  DirectSearch package, which you can freely download from Maple Application Center.

Here is a procedure that does this not only for lists, but also for expressions similar to polynomials. Formal parameters: P - a list or an expression of type `+`,  var is a list of symbols for which the procedure isolates nonlinear members.

IsolateNonlinearTerms:=proc(P::{list,`+`,polynom},var::list)
uses ListTools;
selectremove(p->degree(p)>1 or degree(p)<0, `if`(P::list,FlattenOnce(P),[op(P)]),var);
end proc:


Examples of use:

w:=[[z, y, x, x^3, 1, -5*y], [x*z, x*y, y, 2*x, 1], [x*z, 1/z^2, z, x*y]];
IsolateNonlinearTerms(w, [x,y,z]);
IsolateNonlinearTerms(x^3-3*x*y^2+y^4+x-5*y+4*z-10, [x,y]);
                         

 

 

restart;
pde1 := (y+z)*(diff(u(x, y, z), x))+(z+x)*(diff(u(x, y, z), y))+(x+y)*(diff(u(x, y, z), z)) = 0;
pdsolve(pde1, u(x,y,z), HINT=strip);


See help on  pdsolve  command  for details.


 

restart:  with(plots):

cnsts := [ 4*x1 + x2 <= 12,
             x1 - x2 >= 2,
                  x1 >= 0,
                  x2 >= 0 ];

[4*x1+x2 <= 12, 2 <= x1-x2, 0 <= x1, 0 <= x2]

(1)

feasibleRegion := inequal(cnsts, x1 = 0 .. 4, x2 = -1 .. 2, nolines):
display(feasibleRegion);

 

animate(implicitplot, [c1*x1+2*x2=0,x1=0..4,x2=-1..2, color=red],c1 = -1 .. 0,
                frames = 50, background = feasibleRegion);

 

 


Edit. Below are 2 more options with animation. In the first case, we plot the function of 2 variables for each value of c1. The third coordinate of the highest point gives you the desired maximum. In the second option, we directly plot the maximum as a function of c1:


 

restart:  with(plots):

cnsts := [ 4*x1 + x2 <= 12,
             x1 - x2 >= 2,
                  x1 >= 0,
                  x2 >= 0 ];

[4*x1+x2 <= 12, 2 <= x1-x2, 0 <= x1, 0 <= x2]

(1)

feasibleRegion := inequal(cnsts, x1 = 2 .. 3, x2 = 0 .. 1, nolines):
display(feasibleRegion);
FeasibleRegion:=display(plottools:-transform((x,y)->[x,y,0])(feasibleRegion));
solve({4*x1 + x2 = 12, x1 - x2 = 2});

 

 

{x1 = 14/5, x2 = 4/5}

(2)

animate(plot3d, [c1*x1+2*x2,x1=x2+2..(12-x2)/4,x2=0..4/5, style=surface, color=khaki],c1 = -10 .. 100, frames=90, background=FeasibleRegion);
animate(plot, [[t,'Optimization:-Maximize'(t*x1+2*x2,cnsts)[1],t=-10..c1]],c1 = -10 .. 100, frames=60);

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

Warning, problem appears to be unbounded

 

 

 


 

Download question_on_animate_with_background_new2.mw

Download question_on_animate_with_background_new1.mw


 

restart;
with(plots):
F:=2+r*(b-2)*sqrt(b-1):
A:=implicitplot([b=2,F], b=-1..5, r=-5..10, linestyle=[3,1],color=[black,red], thickness=[0,2], gridrefine=3):
B:=inequal(F<0, b=1..5, r=-5..10, color="LightBlue", optionsexcluded = [color = "Pink"],nolines, numpoints=10000):

T:=textplot([[1.7,7,"F=0"],[2.5,-3,"F=0"],[4,-3,"F<0"],[1.4,9,"F<0"],[3,3,"F>0"]], font=[times,bold,16]):
display(A,B,T, size=[500,500], view=[-1..5,-5..10]);

 

# The red curve  F=0 can also be plotted explicitly 
plot(-2/(b-2)/sqrt(b-1), b=-1..5, -5..10, color=red, thickness=2, discont);

 

 


 

Download plot5.mw

Another way to find the roots of the equation in a certain range is  RootFinding:-Analytic  command. It works numerically, but the use of  identify  command in some cases allows you to get symbolic (i.e. exact) root values. For comparison, the same example was solved by  Student:-Calculus1:-Roots  command (2000 times slower).

 

restart;
Eq:=sin(x)=simplify(convert(sin(Pi/60),radical));

sin(x) = (1/16)*(-2*3^(1/2)+2)*(5+5^(1/2))^(1/2)+(1/16)*2^(1/2)*(5^(1/2)-1)*(1+3^(1/2))

(1)

CodeTools:-Usage(RootFinding:-Analytic(Eq, re=3*Pi/2..5*Pi, im=-1..1));
evalf[12]([%]);
identify(%);

memory used=3.55MiB, alloc change=0 bytes, cpu time=46.00ms, real time=82.00ms, gc time=0ns

 

9.37241808320955, 15.6556033903891, 12.6187304919190, 6.33554518473960

 

[9.37241808321, 15.6556033904, 12.6187304919, 6.33554518474]

 

[(179/60)*Pi, (299/60)*Pi, (241/60)*Pi, (121/60)*Pi]

(2)

ts:=time():
Student:-Calculus1:-Roots(Eq, x=3*Pi/2..5*Pi);
time()-ts;
identify(%%);

Warning, some roots are returned as numeric approximations

 

[6.335545185, 9.372418083, 12.61873049, 15.65560339]

 

187.782

 

[(121/60)*Pi, (179/60)*Pi, (241/60)*Pi, (299/60)*Pi]

(3)

187.782/0.082;

2290.024390

(4)

 


 

Download roots.mw

It follows from the obvious symmetry of your equations that any pair of equal numbers  xa1=xb1   is a solution to this system. In addition, if we plot graphs, then it is clearly visible on them that there are a couple of solutions when the roots are not equal. One of these pairs was indicated by Carl in his reply. Another symmetric solution is found below:


 

restart;

T:=325:

a12:=2305.28444347652 - 9.14490843016421*T + 0.00680052257590234*T^2:

a21:=-6665.24838284836 + 46.0897018087247*T - 0.0694991633494123*T^2:

alfa:=0.3:

x2:=1-x1:

tau12:=a12/T:

tau21:=a21/T:

G12:=exp(-alfa*tau12):

G21:=exp(-alfa*tau21):

lng1:=x2^2*(tau21*(G21/(x1+x2*G21))^2+tau12*(G12/((x2+x1*G12)^2))):

lng2:=x1^2*(tau12*(G12/(x2+x1*G12))^2+tau21*(G21/((x1+x2*G21)^2))):

lnga1:=subs(x1=xa1,lng1):

lngb1:=subs(x1=xb1,lng1):

lnga2:=subs(x1=xa1,lng2):

lngb2:=subs(x1=xb1,lng2):

r1:=lnga1+ln(xa1)=lngb1+ln(xb1);

(1-xa1)^2*(.4966874722/(.4073000470+.5926999530*xa1)^2+.1510891213/(-0.464212743e-1*xa1+1)^2)+ln(xa1) = (1-xb1)^2*(.4966874722/(.4073000470+.5926999530*xb1)^2+.1510891213/(-0.464212743e-1*xb1+1)^2)+ln(xb1)

(1)

r2:=lnga2+ln(1-xa1)=lngb2+ln(1-xb1);

xa1^2*(.1440753718/(-0.464212743e-1*xa1+1)^2+1.219463331/(.4073000470+.5926999530*xa1)^2)+ln(1-xa1) = xb1^2*(.1440753718/(-0.464212743e-1*xb1+1)^2+1.219463331/(.4073000470+.5926999530*xb1)^2)+ln(1-xb1)

(2)

A:=plots:-implicitplot(r1, xa1=0..1, xb1=0..1, color=red, gridrefine=4);
B:=plots:-implicitplot(r2, xa1=0..1, xb1=0..1, color=blue, gridrefine=4);
plots:-display(A,B);

# One of the solutions at the top
fsolve({r1,r2},{xa1=0..0.2,xb1=0.4..0.6});

 

 

 

{xa1 = 0.6850883394e-1, xb1 = .5048762298}

(3)

 


 

Download NRTL_new.mw

 

Maybe a typo has occurred? If we remove the minus in front of  a , and increase the coefficient itself, we get a very similar graph:

restart;
plots:-implicitplot((7.72-7.72*B)*(25.717267500*a) = 662204.4444*B^2, a = 10 .. 50000, B = 0.01 .. 1, color="Blue", thickness=3, tickmarks=[7,10], gridrefine=3, size=[700,400], gridlines, view=[0..60000,0..1]);

             

 

restart;
P:=x->x^4+x^3+a*x^2+sqrt(2)*x+b;
solve(evalc([Re,Im](P(1+I)))=~[0,0]);
solve(eval(P(x),%));

                

 

 

First 14 15 16 17 18 19 20 Last Page 16 of 242