## How to export data from plot...

Hello, everyone!

I have a problem. Colleague send me a figure in maple-file – .mw. How to export data from this plot? I usually use "plottools:-getdata" command, but what strategy should I apply in case of attached file?

Test-example.mw

## How do I approximate a "smooth" piecewise function...

I am currently struggling with a piecewise functions. Seeing as this function is defined over a large number of breakpoints, it would seem that Maple is unable to execute a certain integration (over a particular interval). Since the first derivatives match up at all breakpoints and the piecewise function looks perfectly smooth to my eyes, I figured that I might obtain a sensible result approximating the function by means of a polynomial and integrating based on the resulting approximation. When I try to use minimax from numapprox, however, I just receive an error message to the effect that:

Error, (in evalf/int/CreateProc) function does not evaluate to numeric

I'd very much appreciate it if anyone had any pointers on how I might find a reasonably precise (not necessarily polynomial) approximation for a seemingly smooth piecewise function with a large number of breakpoints.

Many thanks.

Edit: Please find attached my worksheet: worksheet.mw. My goal is to approximate the function Z over the interval 5/3 to 5 by something that is not defined in a piecewise manner.

## What could be wrong with this code?...

I am trying to run this code but this error message is coming up "Error, (in fprintf) number expected for floating point format
". What have I done wrong and how do I correct it.

Thank you and best regards.

```restart;
Digits:=30:

f:=proc(n)
-((sin(x[n]))/(2-sin(x[n])))*(2+sin(x[n]-Pi)):
end proc:

e1:=y[n+1]=h*delta[n]-(1/2)*h^2*(-u^2*sin((1/2)*u)+2*cos((1/2)*u)*u-2*cos(u)*u-4*sin((1/2)*u)+2*sin(u))*f(n)/(u^2*(2*sin((1/2)*u)-sin(u)))+(1/2)*h^2*(u^2*sin((1/2)*u)+2*cos((1/2)*u)*u-4*sin((1/2)*u)+2*sin(u)-2*u)*f(n+1)/(u^2*(2*sin((1/2)*u)-sin(u)))-(1/2)*h^2*(u*sin(u)+2*cos(u)-2)*f(n+1/2)/(u*(2*sin((1/2)*u)-sin(u)))+y[n]:
e2:=h*delta[n+1]=h*delta[n]+h^2*(u*sin((1/2)*u)+cos(u)-1)*f(n)/(u*(2*sin((1/2)*u)-sin(u)))+h^2*(u*sin((1/2)*u)+cos(u)-1)*f(n+1)/(u*(2*sin((1/2)*u)-sin(u)))-h^2*(u*sin(u)+2*cos(u)-2)*f(n+1/2)/(u*(2*sin((1/2)*u)-sin(u))):
e3:=y[n+1/2]=(1/2)*h*delta[n]-(1/8)*h^2*(-u^2*sin((1/2)*u)+4*cos((1/2)*u)*u-4*cos(u)*u-16*sin((1/2)*u)+8*sin(u))*f(n)/(u^2*(2*sin((1/2)*u)-sin(u)))+(1/8)*h^2*(u*sin((1/2)*u)+4*cos((1/2)*u)-4)*f(n+1)/(u*(2*sin((1/2)*u)-sin(u)))-(1/8)*h^2*(u^2*sin(u)+4*cos(u)*u+16*sin((1/2)*u)-8*sin(u)-4*u)*f(n+1/2)/(u^2*(2*sin((1/2)*u)-sin(u)))+y[n]:
e4:=h*delta[n+1/2]=h*delta[n]-(1/2)*h^2*(-u*sin((1/2)*u)+4*cos((1/2)*u)-2*cos(u)-2)*f(n)/(u*(2*sin((1/2)*u)-sin(u)))+(1/2)*h^2*(u*sin((1/2)*u)+4*cos((1/2)*u)-4)*f(n+1)/(u*(2*sin((1/2)*u)-sin(u)))-(1/2)*h^2*(u*sin(u)+2*cos(u)-2)*f(n+1/2)/(u*(2*sin((1/2)*u)-sin(u))):

inx:=0:
ind:=1:
iny:=2:
h:=Pi/4:
n:=0:
omega:=1:
u:=omega*h:
N:=solve(h*p = 8*Pi, p):

c:=1:
for j from 0 to 1 by 1/2 do
t[j]:=inx+j*h:
end do:

vars:=y[n+1/2],y[n+1],delta[n+1/2],delta[n+1]:

step := [seq](eval(x, x=c*h), c=1..N):
printf("%6s%15s%15s%16s%15s%15s%15s\n",
"h","Num.y","Num.z","Ex.y","Ex.z","Error y","Error z");
#eval(<vars>, solve({e||(1..4)},{vars}));

st := time():
for k from 1 to N do

par1:=x[0]=t[0],x[1/2]=t[1/2],x[1]=t[1]:
par2:=y[n]=iny,delta[n]=ind:
res:=eval(<vars>, fsolve(eval({e||(1..4)},[par1,par2]), {vars}));

for i from 1 to 2 do
exy:=eval(2+sin(c*h/2)):
exz:=eval(-sin(c*h/2)+10*epsilon*cos(10*c*h)):
printf("%6.5f%17.9f%15.9f%15.9f%15.9f %13.5g%15.5g\n",
h*c,res[i],res[i+2],exy,exz,abs(res[i]-exy),
abs(res[i+2]-exz)):

c:=c+1:

end do:
iny:=res[2]:
ind:=res[4]:
inx:=t[1]:
for j from 0 to 1 by 1/2 do
t[j]:=inx + j*h:
end do:
end do:```

## How do I plot handle the string in this plot?...

In this code, there are 8 vectors, each vector has 5 points except vectors 7 and 8 with 4 points. I made their  (7 and 8) respective fifth point a string. However, this error statement pops up "Error, invalid input: log[10] expects its 1st argument, x, to be of type algebraic, but received 2.51E-10".

Can someone help me with the correction?

Thank you all and kind regards

restart;

B:=<<601,1201,2401,4801,9601>|<2.87E-19,7.94E-21,7.94E-23,1.03E-24,5.91E-26>|
<2001,4001,8001,16001,32001>|<3.30E-2,5.37E-4,8.47E-6,1.33E-7,2.08E-9>|
<183,365,728,1476,2910>|<4.58E0,7.54E-8,2.20E-11,4.95E-14,9.15E-15>|
<1621,3020,6166,12022,"2222">|<2.95E-3,8.51E-6,3.39E-8,2.51E-10,"2.51E-10">>:

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]):

B[i, 1] := log[10](B[i, 1]):
B[i, 3] := log[10](B[i, 3]):
B[i, 5] := log[10](B[i, 5]):
B[i, 7] := log[10](B[i, 7]):

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","", "BNM","", "BHM","", "ARKN"],
#title="Efficiency Curve for Example 5",
style = ["pointline","point","pointline","point","pointline","point","pointline","point"],
symbolsize = 15,axes = framed,
symbol = [box,box, circle,circle,solidbox, solidbox,solidcircle, solidcircle],
color=[ red, red, gold,gold, blue, blue,black, black],
axis = [gridlines = [colour = green, majorlines = 1,linestyle = dot]],
labels = [typeset(log__10(`NFE`)), typeset(log__10(`Max Err`))]);

## Profiling with sub-procedures...

Hi!

Is there a way to use the profiling function for procedures within a procedure? I have this code as an example:

Test:=proc(a,b)::real;

local c;Ergebnis;

Test2:= proc(d,e) #Prozedur zum Schreiben der Ausgaben

f:=d+e;
g:=f*d;
5+3;
3/0;
end proc:
c:=a+b;
Ergebnis:=Test2(a,c)
end proc:

infolevel[Test]:= 2:
CodeTools:-Profiling:-Profile(Test):
Test(3,4);
CodeTools:-Profiling:-PrintProfiles(Test);

just profiling Test2, the subprocedure, doesn't do anything. Obviously, in this example, the division by zero is the problem, however, I am interested in a general solution to detect more complicated problems within a sub-procedure. Is there a way to also profile it, and thus see where the computation stops?

## Is it possible to tabulate multiple math plots wit...

Hi!

I am trying to construct multiple math plots within one procedure, however, only one is displayed at a time, and it is always at the end of the outputs. I searched for that peculiar behaviour and eventually found on

https://de.maplesoft.com/support/help/Maple/view.aspx?path=DataFrame/Tabulate

the line

This command inserts the assembly of Tables and Embedded Components into the worksheet using the InsertContent facility. The inserted content is placed after any usual output of the Execution Group in which this command is called. Each Execution Group allows for only one inserted result to exist at any given time. Multiple calls to commands which utilize the InsertContent facility made within the same Execution Group will result in each successive inserted assembly replacing any assembly inserted earlier for that Execution Group.

Is there a workaround for this issue? I cannot find one, e.g. on

https://de.maplesoft.com/support/help/Maple/view.aspx?path=DocumentTools%2fInsertContent

## On creating a DataFrame with ArrayTools...

Hi!

I have trouble using ArrayTools, specifically Extend, to create a DataFrame. DF1 and DF2 are, even though awkwardly, created correctly, DF 3 does not work. My code is this:

with(ArrayTools):

with(DocumentTools):

nOben:=10;

nUnten:=3;

InitialisierungDF:=Vector(nOben-nUnten+1); #Initialisierung des Dataframes für die Tabelle
InitialisierungSpalte:=Vector[row](nOben-nUnten+1);
for i from 1 to nOben-nUnten+1 do
InitialisierungDF(i):=oE;
InitialisierungSpalte(i):=n=i;
end do;
print(InitialisierungDF);
print(InitialisierungSpalte);
DF1:= DataFrame( < InitialisierungDF | InitialisierungDF | InitialisierungDF | InitialisierungDF>,
columns = [ GK, PZP, PZM, PY],
rows = InitialisierungSpalte);
DF2:= DataFrame( < InitialisierungDF | InitialisierungDF | InitialisierungDF | InitialisierungDF>,
columns = [ GK, PZP, PZM, PY],
rows = InitialisierungSpalte);
Tabulate(DF1);

print(InitialisierungSpalte);
print(Extend(InitialisierungDF,[oE]));
print(Extend(InitialisierungSpalte,[RMSE]));
InitialisierungDF:=Extend(InitialisierungDF,[oE]);
InitialisierungSpalte:=Extend(InitialisierungSpalte,[RMSE]);
DF3:= DataFrame( < InitialisierungDF | InitialisierungDF | InitialisierungDF | InitialisierungDF>,
columns = [ GK, PZP, PZM, PY],
rows = InitialisierungSpalte);

The result is supposed to be a math table of the form

GK   PZP    PZM   PY

n=1     oE     oE       oE     oE

n=2     oE     oE       oE     oE

...

RMSE  oE    oE      oE       oE

I have run into multiple problems, potentially bugs. First, I have to initialize awkwardly with

InitialisierungDF:=Vector(nOben-nUnten+1); #Initialisierung des Dataframes für die Tabelle
InitialisierungSpalte:=Vector[row](nOben-nUnten+1);
for i from 1 to nOben-nUnten+1 do
InitialisierungDF(i):=oE;
InitialisierungSpalte(i):=n=i;
end do;

because InitialisierungDF:=Vector(1..nOben-nUnten+1,oE) doesn't work and produces entries like oE(1) for whatever reason.

Then, creating the data frame DF3 does not work at all with the above mentioned code.

InitialisierungDF:=Vector(nOben-nUnten+1); #Initialisierung des Dataframes für die Tabelle
InitialisierungSpalte:=Vector[row](nOben-nUnten+1);
for i from 1 to nOben-nUnten+1 do
InitialisierungDF(i):=oE;
InitialisierungSpalte(i):=n=i;
end do;
print(InitialisierungDF);
print(InitialisierungSpalte);
DF1:= DataFrame( < InitialisierungDF | InitialisierungDF | InitialisierungDF | InitialisierungDF>,
columns = [ GK, PZP, PZM, PY],
rows = InitialisierungSpalte);
DF2:= DataFrame( < InitialisierungDF | InitialisierungDF | InitialisierungDF | InitialisierungDF>,
columns = [ GK, PZP, PZM, PY],
rows = InitialisierungSpalte);

print(InitialisierungSpalte);
print(Extend(InitialisierungDF,[oE]));
print(Extend(InitialisierungSpalte,[RMSE]));

DF3:= DataFrame(<Extend(InitialisierungDF,[oE])|Extend(InitialisierungDF,[oE])|Extend(InitialisierungDF),[oE]|Extend(InitialisierungDF,[oE])>,
columns = [ GK, PZP, PZM, PY],
rows = Extend(InitialisierungSpalte,[RMSE]));

doesn't work either, however, this time a different error occurs. Looking at the vectors the dimensions seem to be correct in each case, though. Any suggestions on how to fix these things?

## How to construct a math plot within a procedure?...

Good evening,

is it possible to construct a math table with Maple just with Code? Something like

test_1    test_2

a       9,6        5,6

b     pending  8,4

to be able to print it? There is

https://de.maplesoft.com/support/help/Maple/view.aspx?path=worksheet/documenting/table            and

https://de.maplesoft.com/support/help/Maple/view.aspx?path=worksheet%2freference%2ftableproperties

, however, I do not understand how to just get any table going and assign calculated values within a procedure as is possible with a matrix for example. Thank you in advance! :)

## Code for a math table...

Hi there!

I have a procedure that compares the (2n+1)-point Gauß-Kronrod-Quadrature to the (2n+1)-point Patterson-Quadratures for a range of n. I have plotted the results (the absolute and relative error if they "exist", meaning they need to posess certain features) in a graph, however they do not look very insightful. For the lowest n, the reader gets an impression for the different accuracy of the quadrature rules, however, for higher n's, the resulting points are basically just on the x-axis with no difference to see. Is it possible to also print a math table with Maple? Something like:

GKQ   PZ+   PZ-   PY

1     1,04   1,03  1,02  1,02

2     1,09   1,04     -       -

3     1,02   1,01  1,01  1,01

4     1,03   1,02  1,01  1,01

with - meaning no existance for that particular n? I havent found anything about that on the internet, it's all about plotting.

My (long) code is this:

restart:
with(LinearAlgebra):
with(ListTools):
with(PolynomialTools):
with(CurveFitting):
with(plots):
Plotting:=proc(Unten,Oben,f,g,nUnten,nOben)::plot;

local SpeicherlisteX, SpeichervektorX, #speichert die Stützstellen
SpeichervektorXGekürzt, #streicht nicht existierende Quaraturformeln.
SpeicherlisteYAbs, SpeichervektorYAbs,  #speichert die Stützwerte des späteren Splines aus dem absoluten Quadraturfehler
SpeicherlisteYRel, SpeichervektorYRel,  #speichert die Stützwerte des späteren Splines aus den relativen Quadraturfehler
î, #Laufvariable
InterpolationsfunktionAbs, #speichert den Spline aus dem absoluten Interpolationsfehler
InterpolationsfunktionRel, #speichert den Spline aus den relativen Fehlern von f
GraphAbsGK, GraphAbsPY, GraphAbsPZP, GraphAbsPZM, #speichert den Graphen aus dem Spline aus dem absoluten Interpolationsfehler          GraphRelGK, GraphRelPY, GraphRelPZP, GraphRelPZM, #speichert den Graphen aus dem Spline aus den relativen Fehlern  von f
PunkteAbsGK, PunkteAbsPY, PunkteAbsPZP, PunkteAbsPZM,#speichert den Punktgraphen aus dem absoluten Interpolationsfehler
PunkteRelGK, PunkteRelPY, PunkteRelPZP, PunkteRelPZM, #speichert den Punktgraphen aus dem absoluten Interpolationsfehler
NichtexistenzGK, NichtexistenzPY, NichtexistenzPZP, NichtexistenzPZM, #speichert die Häufigkeit der Nichtexistenz

p,i,c,d,e,Hn,Koeffizienten,s,j,M,V,S,K,nNeu,Em,Hnm,KnotenHnm,KoeffizientenHnm,h0,b,gxi,Gewichte,Delta,Ergebnis,
Endergebnis,Koeffizient,Rest,a,VorgegebeneKnoten,TatsächlicherWert, DoppelterKnoten, KomplexerKnoten,

Text:= proc() #Prozedur zum Schreiben der Ausgabe
uses T= Typesetting;
T:-mrow(seq(`if`(e::string, T:-mn(e), T:-Typeset(T:-EV(e))), e= [args]))
end proc,
OrtPol:= proc(G,N)::list; #Prozedur zum Berechnen der benötigten orthogonalen Polynome
local q,r,R;
q[-1]:=0;
q[0]:=1;

for r from 1 to N do
end do;
return(fsolve(q[N]));
end proc,
BasenwechselNormiert:=proc(Dividend, m)::list; #stellt ein gegebenes Polynom über eine Linearkombination der orthogonalen Polynome #dar.
local BasenwechselNormiert;

Koeffizient:=quo(Dividend, p[m],x);

Rest:=rem(Dividend, p[m],x);

if m=0 then
BasenwechselNormiert:=[Koeffizient*evalf(Int(g*p[m]^2,x=Unten..Oben))];
else

BasenwechselNormiert:=[Koeffizient*evalf(Int(g*p[m]^2,x=Unten..Oben)),op(procname(Rest,m-1))];

end if;

end proc,
Basenwechsel:=proc(Dividend, m)::list; #stellt ein gegebenes Polynom über eine Linearkombination der orthogonalen Polynome 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,
Erweiterung:= proc(Unten, Oben, f,g,Liste,n)::real; #Prozedur zur Berechnung der optimalen Erweiterung nach Knotenvorgabe
#Unten:= Untere Intervallgrenze; Oben:= Obere Intervallgrenze; f:= zu integrierende Funktion;
#g:= Gewicht; Liste:= Liste der alten Knoten, n:= Anzahl hinzuzufügender Knoten;

Hn:=mul(x-Liste[i],i=1..numelems(Liste));

Koeffizienten:=FromCoefficientList(BasenwechselNormiert(Hn,numelems(Liste)+1),x,termorder=reverse); #Die Koeffizienten der orthogonalen Polynome werden hier als Koeffizienten der Monome gespeichert.

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
if s<>j then
M(j+1,s+1):=M(s+1,j+1);
end if;
end do;

end do;

S:=LinearSolve(M,V);
K:=evalindets(S,name,()->2);

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

if (KnotenHnm[1]<-1-10^(-10)) or (KnotenHnm[n+numelems(Liste)]>1+10^(-10)) then
return(false)
else
KomplexerKnoten:=false;
for i from 1 to n+numelems(Liste) do

if(Im(KnotenHnm[i])>10^(-10)) then
KomplexerKnoten:=true
end if;
end do;
if KomplexerKnoten=true then
return(false)
else
DoppelterKnoten:=false;
for i from 1 to n+numelems(Liste)-1 do

if (KnotenHnm[i+1]-KnotenHnm[i]<10^(-10)) then
DoppelterKnoten:=true
end if;
end do;
if DoppelterKnoten=true then
return(false)
else

KoeffizientenHnm:=Reverse(Basenwechsel(Hnm,n+numelems(Liste)));  #Das Polynom Hnm wird über die orthogonalen Polynome dargestellt.

h0:=evalf(Int(g,x=Unten..Oben)); #Beginn der Berechnung der Gewichte

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

b[j]:=KoeffizientenHnm[j+1]+(d[j]+KnotenHnm[i]*c[j])*b[j+1]+e[j+1]*b[j+2];

end do;

gxi:=quo(Hnm,x-KnotenHnm[i],x);

Gewichte[i]:=c[0]*b[1]*h0/eval(gxi,x=KnotenHnm[i]);

Delta[i]:=c[0]*b[1];
end do;

Endergebnis:=Re(evalf(Ergebnis))
end if;
end if;
end if;
end proc:

p[-1]:=0;
p[0]:=1;
for i from 1 to (2*nOben+1)*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

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;

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 2*nOben+1 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 2*nOben+1 do
for j from 2 to s do

end do;
end do;
for î from nUnten to nOben do
VorgegebeneKnoten[î]:=OrtPol(g,î);
end do;
TatsächlicherWert:=evalf(Int(f*g,x= Unten..Oben));
GraphAbsGK:=plot([]); PunkteAbsGK:=plot([]); GraphAbsPZP:=plot([]); PunkteAbsPZP:=plot([]); GraphAbsPZM:=plot([]); PunkteAbsPZM:=plot([]); GraphAbsPY:=plot([]); PunkteAbsPY:=plot([]);
GraphRelGK:=plot([]); PunkteRelGK:=plot([]); GraphRelPZP:=plot([]); PunkteRelPZP:=plot([]); GraphRelPZM:=plot([]); PunkteRelPZM:=plot([]); GraphRelPY:=plot([]); PunkteRelPY:=plot([]);
SpeicherlisteX:=[];
SpeicherlisteYAbs:=[];
SpeicherlisteYRel:=[];
for î from nUnten to nOben do
if Erweiterung(Unten,Oben,f,g,[VorgegebeneKnoten[î]],î+1) <> false then
SpeicherlisteX:=[op(SpeicherlisteX),î]; #Stützstellen definieren
SpeicherlisteYAbs:=[op(SpeicherlisteYAbs),Erweiterung(Unten,Oben,f,g,[VorgegebeneKnoten[î]],î+1)-evalf(Int(f*g,   x=Unten..Oben))]; #Bestimmen des absoluten Fehlers von f für n=î
if abs(TatsächlicherWert) > 10^(-10) then #Bestimmen des relativen Fehlers von f1 falls                                                          #dieser definiert ist
SpeicherlisteYRel:=[op(SpeicherlisteYRel),abs(SpeicherlisteYAbs[-1]/TatsächlicherWert)];
end if;
end if;
end do;
if numelems(SpeicherlisteX)>0 then
SpeichervektorX:=Vector[row](numelems(SpeicherlisteX),SpeicherlisteX);
SpeichervektorYAbs:=Vector[row](numelems(SpeicherlisteYAbs),SpeicherlisteYAbs);
PunkteAbsGK:= plot(SpeichervektorX,SpeichervektorYAbs,style = point, color=red, legend = ["GK"]);
#  Generierung des Punktgraphen, der sich aus den absoluten Fehlern von f ergibt
if numelems(SpeicherlisteX)>1 then
InterpolationsfunktionAbs:=Spline(SpeichervektorX,SpeichervektorYAbs,n);
#  Splines aus Stützpunkten, die sich aus den absoluten Fehlern von f ergeben
GraphAbsGK:= plot(InterpolationsfunktionAbs, n=nUnten..nOben, color=red);
#  Generierung des Graphen, der sich aus dem Spline aus den absoluten Fehlern von f ergibt
end if;
end if;

if abs(TatsächlicherWert) > 10^(-10) then
# falls der relative Fehler definiert ist analoges Vorgehen für die relativen Fehler
if numelems(SpeicherlisteX)>0 then
SpeichervektorYRel:=Vector[row](numelems(SpeicherlisteYRel),SpeicherlisteYRel);
PunkteRelGK:= plot(SpeichervektorX,SpeichervektorYRel,style = point, color=red, legend = ["GK"]);

if numelems (SpeicherlisteX)>1 then
InterpolationsfunktionRel:=Spline(SpeichervektorX,SpeichervektorYRel,n);
GraphRelGK:= plot(InterpolationsfunktionRel, n=nUnten..nOben, color=red);
end if;
end if;
end if;
NichtexistenzGK:=nOben-nUnten+1-numelems(SpeicherlisteX);

SpeicherlisteX:=[]; # analoges Vorgehen für PZP
SpeicherlisteYAbs:=[];
SpeicherlisteYRel:=[];
for î from nUnten to nOben do
if Erweiterung(Unten,Oben,f,g,[-1,VorgegebeneKnoten[î]],î) <> false then
SpeicherlisteX:=[op(SpeicherlisteX),î]; #Stützstellen definieren
SpeicherlisteYAbs:=[op(SpeicherlisteYAbs),Erweiterung(Unten,Oben,f,g,[-1,VorgegebeneKnoten[î]],î)-TatsächlicherWert]; #Bestimmen des absoluten Fehlers von f für n=î
if abs(TatsächlicherWert) > 10^(-10) then #Bestimmen des relativen Fehlers von f1 falls                                                          #dieser definiert ist
SpeicherlisteYRel:=[op(SpeicherlisteYRel),abs(SpeicherlisteYAbs[-1]/TatsächlicherWert)];
end if;
end if;
end do;
if numelems(SpeicherlisteX)>0 then
SpeichervektorX:=Vector[row](numelems(SpeicherlisteX),SpeicherlisteX);
SpeichervektorYAbs:=Vector[row](numelems(SpeicherlisteYAbs),SpeicherlisteYAbs);
PunkteAbsPZP:= plot(SpeichervektorX,SpeichervektorYAbs,style = point, color=orange);
#  Generierung des Punktgraphen, der sich aus den absoluten Fehlern von f ergibt
if numelems(SpeicherlisteX)>1 then
InterpolationsfunktionAbs:=Spline(SpeichervektorX,SpeichervektorYAbs,n);
#  Splines aus Stützpunkten, die sich aus den absoluten Fehlern von f ergeben
GraphAbsPZP:= plot(InterpolationsfunktionAbs, n=nUnten..nOben, color=orange);
#  Generierung des Graphen, der sich aus dem Spline aus den absoluten Fehlern von f ergibt
end if;
end if;

if abs(TatsächlicherWert) > 10^(-10) then
# falls der relative Fehler definiert ist analoges Vorgehen für die relativen Fehler
if numelems(SpeicherlisteX)>0 then
SpeichervektorYRel:=Vector[row](numelems(SpeicherlisteYRel),SpeicherlisteYRel);
PunkteRelPZP:= plot(SpeichervektorX,SpeichervektorYRel,style = point, color=orange, legend = ["PZP"]);

if numelems (SpeicherlisteX)>1 then
InterpolationsfunktionRel:=Spline(SpeichervektorX,SpeichervektorYRel,n);
GraphRelPZP:= plot(InterpolationsfunktionRel, n=nUnten..nOben, color=orange);
end if;
end if;
end if;
NichtexistenzPZP:=nOben-nUnten+1-numelems(SpeicherlisteX);
SpeicherlisteX:=[];# analoges Vorgehen für PZM
SpeicherlisteYAbs:=[];
SpeicherlisteYRel:=[];
for î from nUnten to nOben do
if Erweiterung(Unten,Oben,f,g,[VorgegebeneKnoten[î],1],î) <> false then
SpeicherlisteX:=[op(SpeicherlisteX),î]; #Stützstellen definieren
SpeicherlisteYAbs:=[op(SpeicherlisteYAbs),Erweiterung(Unten,Oben,f,g,[VorgegebeneKnoten[î],1],î)-TatsächlicherWert]; #Bestimmen des absoluten Fehlers von f für n=î
if abs(TatsächlicherWert) > 10^(-10) then #Bestimmen des relativen Fehlers von f1 falls                                                          #dieser definiert ist
SpeicherlisteYRel:=[op(SpeicherlisteYRel),abs(SpeicherlisteYAbs[-1]/TatsächlicherWert)];
end if;
end if;
end do;
if numelems(SpeicherlisteX)>0 then
SpeichervektorX:=Vector[row](numelems(SpeicherlisteX),SpeicherlisteX);
SpeichervektorYAbs:=Vector[row](numelems(SpeicherlisteYAbs),SpeicherlisteYAbs);
PunkteAbsPZM:= plot(SpeichervektorX,SpeichervektorYAbs,style = point, color=blue, legend = ["PZM"]);
#  Generierung des Punktgraphen, der sich aus den absoluten Fehlern von f ergibt
if numelems(SpeicherlisteX)>1 then
InterpolationsfunktionAbs:=Spline(SpeichervektorX,SpeichervektorYAbs,n);
#  Splines aus Stützpunkten, die sich aus den absoluten Fehlern von f ergeben
GraphAbsPZM:= plot(InterpolationsfunktionAbs, n=nUnten..nOben, color=blue);
#  Generierung des Graphen, der sich aus dem Spline aus den absoluten Fehlern von f ergibt
end if;
end if;

if abs(TatsächlicherWert) > 10^(-10) then
# falls der relative Fehler definiert ist analoges Vorgehen für die relativen Fehler
if numelems(SpeicherlisteX)>0 then
SpeichervektorYRel:=Vector[row](numelems(SpeicherlisteYRel),SpeicherlisteYRel);
PunkteRelPZM:= plot(SpeichervektorX,SpeichervektorYRel,style = point, color=blue, legend = ["PZM"]);

if numelems (SpeicherlisteX)>1 then
InterpolationsfunktionRel:=Spline(SpeichervektorX,SpeichervektorYRel,n);
GraphRelPZM:= plot(InterpolationsfunktionRel, n=nUnten..nOben, color=blue);
end if;
end if;
end if;
NichtexistenzPZM:=nOben-nUnten+1-numelems(SpeicherlisteX);
SpeicherlisteX:=[]; #analoges Vorgehen für PY
SpeicherlisteYAbs:=[];
SpeicherlisteYRel:=[];
for î from nUnten to nOben do
if Erweiterung(Unten,Oben,f,g,[-1,VorgegebeneKnoten[î],1],î-1) <> false then
SpeicherlisteX:=[op(SpeicherlisteX),î]; #Stützstellen definieren
SpeicherlisteYAbs:=[op(SpeicherlisteYAbs),Erweiterung(Unten,Oben,f,g,[-1,VorgegebeneKnoten[î],1],î-1)-TatsächlicherWert]; #Bestimmen des absoluten Fehlers von f für n=î
if abs(TatsächlicherWert) > 10^(-10) then #Bestimmen des relativen Fehlers von f1 falls                                                          #dieser definiert ist
SpeicherlisteYRel:=[op(SpeicherlisteYRel),abs(SpeicherlisteYAbs[-1]/TatsächlicherWert)];
end if;
end if;
end do;
if numelems(SpeicherlisteX)>0 then
SpeichervektorX:=Vector[row](numelems(SpeicherlisteX),SpeicherlisteX);
SpeichervektorYAbs:=Vector[row](numelems(SpeicherlisteYAbs),SpeicherlisteYAbs);
PunkteAbsPY:= plot(SpeichervektorX,SpeichervektorYAbs,style = point, color=purple, legend = ["PY"]);
#  Generierung des Punktgraphen, der sich aus den absoluten Fehlern von f ergibt
if numelems(SpeicherlisteX)>1 then
InterpolationsfunktionAbs:=Spline(SpeichervektorX,SpeichervektorYAbs,n);
#  Splines aus Stützpunkten, die sich aus den absoluten Fehlern von f ergeben
GraphAbsPY:= plot(InterpolationsfunktionAbs, n=nUnten..nOben, color=purple);
#  Generierung des Graphen, der sich aus dem Spline aus den absoluten Fehlern von f ergibt
end if;
end if;

if abs(TatsächlicherWert) > 10^(-10) then
# falls der relative Fehler definiert ist analoges Vorgehen für die relativen Fehler
if numelems(SpeicherlisteX)>0 then
SpeichervektorYRel:=Vector[row](numelems(SpeicherlisteYRel),SpeicherlisteYRel);
PunkteRelPY:= plot(SpeichervektorX,SpeichervektorYRel,style = point, color=purple, legend = ["PY"]);

if numelems (SpeicherlisteX)>1 then
InterpolationsfunktionRel:=Spline(SpeichervektorX,SpeichervektorYRel,n);
GraphRelPY:= plot(InterpolationsfunktionRel, n=nUnten..nOben, color=purple);
end if;
end if;
end if;
NichtexistenzPY:=nOben-nUnten+1-numelems(SpeicherlisteX);
print(display({GraphAbsGK,PunkteAbsGK,GraphAbsPZP,PunkteAbsPZP, GraphAbsPZM,PunkteAbsPZM, GraphAbsPY,PunkteAbsPY}, title= "Absoluter Fehler", titlefont=["ROMAN",18]));
if abs(TatsächlicherWert) > 10^(-10) then
print(display({GraphRelGK,PunkteRelGK,GraphRelPZP,PunkteRelPZP, GraphRelPZM,PunkteRelPZM, GraphRelPY,PunkteRelPY}, title= "Relativer Fehler", titlefont=["ROMAN",18]));
end if;
Text("Häufigkeit der Nichtexistenz: GK ",NichtexistenzGK, ", PZP ",NichtexistenzPZP, ", PZM ", NichtexistenzPZM, ", PY ", NichtexistenzPY);

end proc

An example of how it should not look like is this:

Plotting(-1,1,2*x^2+2,1,3,10)

On a side note, Maple's return is

Warning, `GraphRelGK` is implicitly declared local to procedure `Plotting`
Warning, `GraphRelPZP` is implicitly declared local to procedure `Plotting`
Warning, `GraphRelPZM` is implicitly declared local to procedure `Plotting`
Warning, `GraphRelPY` is implicitly declared local to procedure `Plotting`

even though I did declare them.

Any suggestions on that minor issue? And on how to construct a math table which allows for a symbol like - for nonexistence?

## How do I solve the attached problem in maple?...

Hi,

I would be really grateful if someone can help me in solving the below attached problem in maple.

## MAPLE is taking too long time...

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

## NewsVendor Problem solution for general Demand dis...

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?

## On seeing what Maple is taking so long for...

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

## How do I make this plot code run?...

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

```

## Boolean algebra in combination with "for i from 1 ...

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