Maple 2016 Questions and Posts

These are Posts and Questions associated with the product, Maple 2016

Hi,

I was trying to differentiate an equation with two random variables with respect to my decision variable. But it took 4-5 hrs still MAPLE could not evaluate nor show an error.

 Later I want to find q* from the FOC and wish to substitute it back into the equation. I will be grateful if someone please help me in doing both.

doubt_1.mw

Thanks in Advance

Hi, 
I want to solve a simple newsvendor problem for general Demand Probability Distribution with maple code.
let, 
Underage cost 
                              c[u]

Overage cost 
                              c[o]

Demand (
D;
) follows pdf 
f(.);
 and CDF 
F;

(.);


Decision varible is quantity 
q;


So cost function is 

TC := piecewise(q-D >= 0, c[o]*(q-D), q-D < 0, c[u]*(-q+D));
 piecewise(0 <= q - D, c[o] (q - D), q - D < 0, c[u] (-q + D))

After taking expectation of TC and derivating it w.r.t q we will get 
q;
* as critical fractile i.e.

q*= F^(-1) (
                             c[u]    
                          -----------
                          c[u] + c[o]
)      i.e. F inverse (
c[u]/(c[u]+c[o]);
)

Can someone please help me how to find q* in MAPLE using solve() or any other function?

 

Thanks in Advance!

Hi!

I have the following code to calculate the optimal quadrature rule with preassigned nodes (within a list), while the number of nodes that are to be added is n. The calculated quadrature rule is then used to approximate an integral.

 

restart;
with(LinearAlgebra):     
with(ListTools):
with(PolynomialTools):
Erweiterung:= proc(Unten, Oben, f,g,Liste,n)::real;
  #Unten:= Untere Intervallgrenze; Oben:= Obere Intervallgrenze; f:= zu integrierende Funktion;
  #G:= Gewicht; Liste:= Liste der alten Knoten, n:= Anzahl hinzuzufügender Knoten;
  local Basenwechsel,p,i,c,d,e,Hn,Koeffizienten,a,s,j,M,V,S,K,nNeu,Em,Hnm,KnotenHnm,KoeffizientenHnm,h0,b,gxi,Gewichte,Delta,Ergebnis,Endergebnis,Koeffizient,Rest;
  uses LinearAlgebra, ListTools;
  Basenwechsel1:=proc(Dividend, m)::list; #stellt ein gegebenes Polynom über eine Linearkombination der orthogonalen Polynomen dar.
   local Basenwechsel;
 
  Koeffizient:=quo(Dividend, p[m],x);

  Rest:=rem(Dividend, p[m],x);
 
  if m=0 then
    Basenwechsel:=[Koeffizient*evalf(int(g*p[m]^2,x=Unten..Oben))];
  else

    Basenwechsel:=[Koeffizient*evalf(int(g*p[m]^2,x=Unten..Oben)),op(procname(Rest,m-1))];
   
  end if;
 
  end proc;
  Basenwechsel2:=proc(Dividend, m)::list; #stellt ein gegebenes Polynom über eine Linearkombination der orthogonalen Polynomen dar.
   local Basenwechsel;
 
  Koeffizient:=quo(Dividend, p[m],x);

  Rest:=rem(Dividend, p[m],x);
 
  if m=0 then
    Basenwechsel:=[Koeffizient];
  else

    Basenwechsel:=[Koeffizient,op(procname(Rest,m-1))];
   
  end if;
 
  end proc;
p[-1]:=0;
p[0]:=1;
for i from 1 to (numelems(Liste)+n)*2 do
  p[i]:=(x^i-add(evalf(int(x^i*p[j]*g,x=Unten..Oben))*p[j]/evalf(int(p[j]^2*g,x=Unten..Oben)),j=0..i-1)); #Berechnung einer Folge orthogonaler Polynome bezüglich der gegebenen Gewichtsfunktion und des gegebenes Intervalles
  pprint(p[i]);
c[i-1]:=coeff(p[i],x,i)/coeff(p[i-1],x,i-1); #Berechnung der dreigliedrigen Rekursion der errechneten orthogonalen Polynome
d[i-1]:=(coeff(p[i],x,(i-1))-c[i-1]*coeff(p[i-1],x,(i-2)))/coeff(p[i-1],x,(i-1));
if i <> 1 then
  e[i-1]:=coeff(p[i]-(c[i-1]*x+d[i-1])*p[i-1],x,i-2)/coeff(p[i-2],x,i-2);
else
  e[i-1]:=0;
end if;
end do;
pprint(Liste[1],numelems(Liste));
Hn:=mul(x-Liste[i],i=1..numelems(Liste));
pprint(Hn);
 Koeffizienten:=FromCoefficientList(Basenwechsel1(Hn,numelems(Liste)+1),x,termorder=reverse); #Die Koeffizienten der orthogonalen Polynome werden hier als Koeffizienten der Monome gespeichert.
pprint(Koeffizienten,HIER);

pprint(c,d,e);
a[0][0]:=1; #Beginn der Berechnung der orthogonalen Produkterweiterungen, die Koeffizienten der orthogonalen Polynome werden wieder über die Monome gespeichert (2*x^2+2 bedeutet bspw. [2,0,2,0,0...] für die Koeffizienten)
a[1][0]:=x;
a[1][1]:=-e[1]*c[0]/c[1]+(d[0]-d[1]*c[0]/c[1])*x+c[0]/c[1]*x^2;
for s from 2 to numelems(Liste)+n do
  a[s][0]:=x^s;
  a[s][1]:=-e[s]*c[0]/c[s]*x^(s-1)+(d[0]-d[s]*c[0]/c[s])*x^s+c[0]/c[s]*x^(s+1);
    pprint (coeff(a[s][1],x,s),as1s);
end do;
for s from 2 to numelems(Liste)+n do
  for j from 2 to s do
    
      pprint(c[j-1]*sum(coeff(a[s][j-1],x,k-1)/c[k-1]*x^k,k=abs(s-j)+2..s+j));  pprint(sum((d[j-1]-c[j-1]*d[k]/c[k])*coeff(a[s][j-1],x,k)*x^k,k=abs(s-j)+1..s+j-1));  pprint(c[j-1]*sum(e[k+1]*coeff(a[s][j-1],x,k+1)/c[k+1]*x^k,k=abs(s-j)..s+j-2));pprint(e[j-1]*sum(coeff(a[s][j-2],x,k)*x^k,k=s-j+2..s+j-2));

     a[s][j]:=c[j-1]*sum(coeff(a[s][j-1],x,k-1)/c[k-1]*x^k,k=abs(s-j)+2..s+j)+sum((d[j-1]-c[j-1]*d[k]/c[k])*coeff(a[s][j-1],x,k)*x^k,k=abs(s-j)+1..s+j-1)-c[j-1]*sum(e[k+1]*coeff(a[s][j-1],x,k+1)/c[k+1]*x^k,k=abs(s-j)..s+j-2)+e[j-1]*sum(coeff(a[s][j-2],x,k)*x^k,k=abs(s-j)+2..s+j-2);

      
   
    
  end do;
end do;
for s from 0 to numelems(Liste)+n do
  for j from 0 to s do
    pprint(a[s][j], Polynom[s][j]);
  end do;
end do;
M:=Matrix(n,n); #Beginn der Erstellung eines linearen Gleichungssystems, dessen Lösung die Koeffizienten der orthogonalen Polynome sind, deren Summe Em die hinzuzufügenden Knoten als Nullstellen hat.
V:=Vector(n);
 
  for s from 0 to n-1 do
    for j from 0 to s do
      M(s+1,j+1):=sum(coeff(a[s][j],x,k)*coeff(Koeffizienten,x,k),k=0..n);
      if s<>j then
        M(j+1,s+1):=M(s+1,j+1);
      end if;
      pprint(M,1);
    end do;
    pprint(testb1);pprint(coeff(a[n][s],x,2));pprint(coeff(Koeffizienten,x,2));
    pprint(testb2); pprint(Koeffizienten);
    M(s+1,n+1):=sum(coeff(a[n][s],x,k)*coeff(Koeffizienten,x,k),k=0..n);
    
    pprint(M,V);
  end do;
pprint(M,V);
S:=LinearSolve(M,V);
K:=evalindets(S,name,()->2);
pprint(K,LinSolve);

Em:=add(p[i]*K[i+1],i=0..n); #Erstellen von Em, dessen Nullstellen die hinzuzufügenden Knoten sind
Hnm:=Hn*Em; #Erstellen von Hnm, welches alle Knoten als Nullstelle besitzt
KnotenHnm:=fsolve(Hnm,complex); #Knotenberechnung

 


   
pprint(Hn,alt,Em,neu,Hnm);
pprint(Testergebnis,nNeu);
pprint(fsolve(Hnm),fsolve(nNeu));
KoeffizientenHnm:=Reverse(Basenwechsel2(Hnm,n+numelems(Liste)));  #Das Polynom Hnm wird über die orthogonalen Polynome dargestellt.
pprint(KoeffizientenHnm);
h0:=int(g,x=Unten..Oben); #Beginn der Berechnung der Gewichte
 pprint(h0,HO);
b[n+numelems(Liste)+2]:=0;
b[n+numelems(Liste)+1]:=0;
  for i from 1 to nops([KnotenHnm]) do
    for j from n+numelems(Liste) by -1 to 1 do
      pprint(test21,KnotenHnm,Hnm);
      b[j]:=KoeffizientenHnm[j+1]+(d[j]+KnotenHnm[i]*c[j])*b[j+1]+e[j+1]*b[j+2];
      pprint(b[j]);
    end do;
    pprint(test23);
    gxi:=quo(Hnm,x-KnotenHnm[i],x);
   pprint(gxi);
    Gewichte[i]:=c[0]*b[1]*h0/eval(gxi,x=KnotenHnm[i]);
    pprint(b);
   
    Delta[i]:=c[0]*b[1];
  end do;
print(DieKnoten,KnotenHnm);
print(DieGewichte, Gewichte);
Ergebnis:=add(eval(f,x=KnotenHnm[k])*Gewichte[k],k=1..nops([KnotenHnm]));
Endergebnis:=evalf(Ergebnis)
end proc


The problem is that the code takes very, very long to run if the weight function is not a polynomial.

Erweiterung(-1,1,2*x^2+2,1,[-0.906179845938664],4)

for example, is done immediately (1, the 4th entry, being the weight function), while

Erweiterung(-1,1,2*x^2+2,2/sqrt(1-x^2),[-.8660254037, 0, .8660254037],4)


takes ages to finish. Is there a tool for me to see what exactly is taking Maple so long? Is there an easy fix, such as evalf()'ing key calculations (other than using (2*x^2+2)*2/sqrt(1-x^2) as the integrand and 1 as the weight function, since the quadrature rules I am looking at are supposed to be good with certain weird weight functions)? Thank you in advance!

Please I have an issue with the attached plot code. Can you kindly help to correct it? 

restart:
interface(rtablesize=infinity):
B:=<<0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.49,"0","0.05","0.1","0.15","0.2","0.25","0.3","0.35","0.4","0.45","0.49">|
	<14.73,14.4,14,13.4,12.67,11.67,10.4,8.67,6,3,0,"14.73","14.4","14","13.4","12.67","11.67","10.4","8.67","6","3","0">|
     <-0.007072,0.013309,0.033707,0.054125,0.074571,0.095056,0.115597,0.136218,0.156956,0.177867,0.199036,0.22059,0.242719,0.265702,0.289932,0.31592,0.344214,0.375124,0.408175,0.441761,0.473484,0.501857>|
	<1.34E+01,1.33E+01,1.33E+01,1.32E+01,1.31E+01,1.31E+01,1.30E+01,1.29E+01,1.28E+01,1.26E+01,1.25E+01,1.22E+01,1.19E+01,1.14E+01,1.08E+01,9.86E+00,8.58E+00,6.90E+00,4.90E+00,2.81E+00,1.00E+00,-2.86E-01>>:

B:
 
 plot([B[..,[1, 2]],B[1..1,[1, 2]], B[.., [3, 4]],B[1..1,[3, 4]], B[..,[5, 6]],B[1..1,[5, 6]],B[.., [7, 8]],B[1..1,[7, 8]],
 	  B[..,[9, 10]],B[1..1,[9, 10]], B[.., [11, 12]],B[1..1,[11, 12]],B[..,[13, 14]],B[1..1,[13, 14]],B[.., [15, 16]],B[1..1,[15, 16]],
 	  B[..,[17, 18]],B[1..1,[17, 18]], B[.., [19, 20]],B[1..1,[19, 20]],B[..,[21, 22]],B[1..1,[21, 22]]],
 	  legend = ["","Experimental","","Simulation"],
 	  style = ["line","line","line","line","line","line","line","line","line","line","line","line","line","line","line","line",
 	 		"line","line","line","line","line","line"],
 	  color=[blue,red], labels=[`V (V)`, `Jsc (mA/cm^2)`]);
 



Download Graph_Example.mw

Hi there!

I am trying to tell Maple to print something one time once a condition has been met for a variable that ranges from a to b, let's say from 1 to 9 in this case. Testbestanden is to be printed, if (F[i]+F[i+1]<20) is true for every i, and Durchgefallen is to be printed as soon as one counterexample is found. The following code illustrates what I need, it doesn't work, however. What is the correct syntax? Switching the lines "if (F[i]+F[i+1]<20)" and "for i from 1 to 9" leads to Maple printing "Testbestanden" for every i that fulfils the statement which is not what I need either.

for i from 1 to 10 do
  F[i]:=i
end do;

if (F[i]+F[i+1]<20)
for i from 1 to 9
then print(Testbestanden)

else print(Durchgefallen)
end if
end do

Thank you in advance! :)

Hello
I solved an equation using Maple . but in answer there is RootOf statement that I can't change to a simple statement without RootOf.

restart

``

``

((Lr2*(vo+vs)*sqrt(Lr1*Lr2)+vo*Lr1*sqrt(Lr1*Lr2))*sin(t*sqrt((Lr1+Lr2)/(Lr1*Lr2*cr)))*sqrt(cr*(Lr1+Lr2))+vs*t*Lr1*(Lr1+Lr2))/((Lr1+Lr2)^2*Lr1)

((Lr2*(vo+vs)*(Lr1*Lr2)^(1/2)+vo*Lr1*(Lr1*Lr2)^(1/2))*sin(t*((Lr1+Lr2)/(Lr1*Lr2*cr))^(1/2))*(cr*(Lr1+Lr2))^(1/2)+vs*t*Lr1*(Lr1+Lr2))/((Lr1+Lr2)^2*Lr1)``

(1)

eq := ((Lr2*(vo+vs)*sqrt(Lr1*Lr2)+vo*Lr1*sqrt(Lr1*Lr2))*sin(t*sqrt((Lr1+Lr2)/(Lr1*Lr2*cr)))*sqrt(cr*(Lr1+Lr2))+vs*t*Lr1*(Lr1+Lr2))/((Lr1+Lr2)^2*Lr1)

((Lr2*(vo+vs)*(Lr1*Lr2)^(1/2)+vo*Lr1*(Lr1*Lr2)^(1/2))*sin(t*((Lr1+Lr2)/(Lr1*Lr2*cr))^(1/2))*(cr*(Lr1+Lr2))^(1/2)+vs*t*Lr1*(Lr1+Lr2))/((Lr1+Lr2)^2*Lr1)

(2)

solve(eq, t)

-(Lr1*Lr2)^(1/2)*(Lr1*cr+Lr2*cr)^(1/2)*sin(RootOf(Lr1*sin(_Z)*(Lr1*cr+Lr2*cr)^(1/2)*((Lr1+Lr2)/(Lr1*Lr2*cr))^(1/2)*(Lr1*Lr2)^(1/2)*vo+Lr2*sin(_Z)*(Lr1*cr+Lr2*cr)^(1/2)*((Lr1+Lr2)/(Lr1*Lr2*cr))^(1/2)*(Lr1*Lr2)^(1/2)*vo+Lr2*sin(_Z)*(Lr1*cr+Lr2*cr)^(1/2)*((Lr1+Lr2)/(Lr1*Lr2*cr))^(1/2)*(Lr1*Lr2)^(1/2)*vs+_Z*Lr1^2*vs+_Z*Lr1*vs*Lr2))*(Lr1*vo+Lr2*vo+Lr2*vs)/(Lr1*(Lr1+Lr2)*vs)

(3)

``

Download spr.mw

How can I simplify to an answer without RootOf?

Thank u

Hi
I want to say maple convert this radicals to a single radical expr.

cos(sqrt(Lr1+Lr2)*t/(sqrt(Cr)*sqrt(Lr1)*sqrt(Lr2)))

convert to :

cos(sqrt((Lr1+Lr2)/(Cr*Lr1*Lr2))*t)
thank u .

In the following code, I want the base of the logarithm in the label to appear at the conventional position not beside the logarithm. How can I achieve this? Thank you all in anticipation of your educative response and do stay safe.

restart;

B:=<<39,76,151,301,601>|<7.71E-8,5.43E-9,3.55E-10,2.11E-11,1.32E-12>|
	<26,51,101,201,401>|<6.46E-3,1.17E-4,1.88E-6,2.96E-8,4.46E-10>|
	<26,51,101,201,401>|<2.74E-4,6.34E-6,1.16E-7,1.85E-9,2.92E-11>|
	<26,51,101,201,401>|<6.48E-4,4.39E-5,2.99E-6,1.88E-7,1.18E-8>>:

for i from 1 to 5 do
   B[i, 2] := log[10](B[i, 2]):		  			
   B[i, 4] := log[10](B[i, 4]): 		
   B[i, 6] := log[10](B[i, 6]): 		
   B[i, 8] := log[10](B[i, 8]):	
       

   
end do:  # computing the log of the max-error
B: # This is the table of values we'll plot.

T:=plot([B[..,[1, 2]],B[1..1,[1, 2]], B[.., [3, 4]],B[1..1,[3, 4]], 
	 B[..,[5, 6]],B[1..1,[5, 6]],B[.., [7, 8]],B[1..1,[7, 8]]], 
	legend = ["","BFFM","", "BHT","", "BHTRKNM", "", "BNM"],
	#title="Efficiency Curve for Example 1",
	style = ["pointline","point","pointline","point","pointline","point","pointline","point"], 
	symbolsize = 15,axes = framed, 
	symbol = [box,box, circle,circle,diamond,diamond,solidcircle, solidcircle],
	color=[red, red, blue,blue, black, black, green, green], 
	axis = [gridlines = [colour = green, majorlines = 1,linestyle = dot]], 
	labels = ["NFE", "log[10](Max Err)"]);

 

I have written the following attached code to use Euler explicit method to solve the following IVP

diff(y(x), x) = 2*(1+x)-y(x), y(2) = 5
With Exact solution  y(x) = 2*x+exp(-x)/exp(-2)

However, I found out that my exact results are not correct while the numerical results are okay. What have I done wrong in the code? Can someone modify the code?

Thank you and kind regards.

Hi all, we know Maple provided discrim method to find discriminant of a polynomial 

I want to write a similar method with independent variable is ,... my code is below

Some examples

Maple already similar method? If not have, we can improve performance it?

Thank you very much.

 

Hi,

I have a random variable that follows Uniform(1,4). Now I have a function which is of the following type:

g := a*alpha+b*t/alpha+exp(alpha)

where,

A := RandomVariable(Uniform(c, d));
                 RandomVariable(Uniform(c, d))
f := proc (alpha) options operator, arrow; PDF(A, alpha) end proc;
alpha -> PDF(A, alpha)
#Defining expectation fuction
E := proc (alpha) options operator, arrow; int(alpha*f(alpha), alpha = c .. d) end proc;
alpha -> int(alpha f(alpha), alpha = c .. d)
#g is a function of random variable &alpha;, where a and b are parameter

 

now I want to find the expectation of g and the first derivative of expectation of g,

E(g)

diff(E(g), t)

 

I understand the way I have defined E(alpha) is improper. But please understand my intent and help! here is the maple file also doubt_1.mw

Dear all,

Please I want only 8 points to show on this curve, how do I specify it?

plot(ln(1+sin(Pi*x)), x = 0 .. 1, legend = numerical, style = point, symbol = box, color = blue, symbolsize = 15, numpoints = 8);

Thank you all and kind regards.

Please do keep safe amidst this global pandemic.

In my worksheet today my intention was to compare the least squares linear regression for three datasets as indicated, but when I right click on the output as seen in the bottom line to select the plot type, all options state there to be independant variables K[0] and K[1], where as the output displays only the variable K as I intended, which part of my code is creating this confusion for maple?

 

 

 

Worksheet Specific Investigation Content

 

S[0] := proc (N, K) options operator, arrow; map(simplify, {seq(seq(seq(piecewise((a^`&varphi;`(b))^(1/(c+1))-floor((a^`&varphi;`(b))^(1/(c+1))) = 0, [a, b, c], NULL), a = 1 .. N), b = 1 .. N), c = 1 .. K)}, 'radical') end proc

T := proc (N, K) options operator, arrow; {seq(seq(seq([a, b, c], a = 1 .. N), b = 1 .. N), c = 1 .. K)} end proc:

S[1] := proc (N, K) options operator, arrow; `minus`(T(N, K), S[0](N, K)) end proc:

CardRatio := proc (N, K) options operator, arrow; nops(S[0](N, K))/nops(S[1](N, K)) end proc:

{CurveFitting[LeastSquares]([seq([k, CardRatio(2, k)], k = 1 .. 10)], K), CurveFitting[LeastSquares]([seq([k, CardRatio(3, k)], k = 1 .. 10)], K), CurveFitting[LeastSquares]([seq([k, CardRatio(4, k)], k = 1 .. 10)], K)}

{1, 44268857/45401356-(532409481/9988298320)*K, 24308311919/13309971675-(135902619982/773879781675)*K}

(1.1)

``

 

 

 

 

Download ask_maple.mw

 

 

Hi,

I am writing the following code and MAPLE is giving me operator error. Please help (file: doubt_6.mw)

d2:=100000000

for m in set_m do
    for n in set_n do
        SOL1 := fsolve({ODE11, ODE12}, {N, t__2});
        N1:=eval(N,SOL1);
        t_2_1 :=eval(t__2,SOL1);
        T_1:= eval(T, [lambda = 3, a = 300, b = .15, c = .25, A__m = 300, A__d = 150, A__r = 50,C__m = 4, P__m = 8, P__d = 10, 
        P__r = 12, theta__m = .15, theta__d = .12, theta__r = 0.5e-1, h__m = .2, h__d = .3, h__r = .5, i__m = .1, i__d = .1, 
        i__r = .1, i__om = .1, i__OD = .15, i__c = .3, i__e = .2, M = 2, alpha = 0.2e-1, t__2 = t_2_1]):
        t_31:= T_1 /m ;
        t_41:= T_1 /(m*n) ;
        if (N1<=t_41 and 2>=t_31) then
            d1:= eval(TCS__1, [lambda = 3, a = 300, b = .15, c = .25, A__m = 300, A__d = 150, A__r = 50, C__m = 4, P__m = 8, 
            P__d = 10, P__r = 12, theta__m = .15, theta__d = .12, theta__r = 0.5e-1, h__m = .2, h__d = .3, h__r = .5, i__m = .1, 
            i__d = .1, i__r = .1, i__om = .1, i__OD = .15, i__c = .3, i__e = .2, M = 2, alpha = 0.2e-1, t__2=t_2_1, N=N1]):
            if (d1<= d2) then
                d2:= d1;
                print("Value is updated",d2,N1,t_2_1,"for",m,n)
            end if
        end if    
        print(N1,t_2_1,t_31,t_41,d1,m,n)
    end do
end do

 

Thanks in advance

 

[[p__jb = (Typesetting[delayDotProduct](c__a . ((rho*(-1+alpha)*t__a-alpha*t__b)/(t__b*t__a)), t__b, true)*t__a+((c__b+`p__-jb`)*alpha-c__b-t__b-`p__-jb`)*t__a-alpha*t__b*(c__b+`p__-jb`))/((2*alpha-2)*t__a-2*alpha*t__b)]]

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