Good day, all.

Please I want to solve the following delay differential equation:

ODE := diff(y(t), t$2) = (2*(1-y(t-1)^2))*(diff(y(t), t))-y(t)

ics := y(0) = 1, (D(y))(0) = 0

using the following codes but there is an error. Please kindly help to modify the codes.

restart;
Digits:=30:

f:=proc(n)
	2*(1-(y[n-1])^2)*delta[n]+y[n]:
end proc:

g:=proc(n)
	-4*y[n-1]*delta[n-1]+2*(1-(y[n-1])^2)*f(n)-delta[n]:
end proc:


e1:=y[n+2] = -y[n]+2*y[n+1]+(1/120)*h^2*(-3*h*g(n+2)+3*g(n)*h+16*f(n+2)+16*f(n)+88*f(n+1)):
e2:=h*delta[n] = -y[n]+y[n+1]-(1/1680)*h^2*(-128*h*g(n+1)-11*h*g(n+2)+59*g(n)*h+40*f(n+2)+520*f(n)+280*f(n+1)):
e3:=h*delta[n+1] = -y[n]+y[n+1]+(1/1680)*h^2*(-152*h*g(n+1)-10*h*g(n+2)+32*g(n)*h+37*f(n+2)+187*f(n)+616*f(n+1)):
e4:=h*delta[n+2] = -y[n]+y[n+1]+(1/1680)*h^2*(128*h*g(n+1)-101*h*g(n+2)+53*g(n)*h+744*f(n+2)+264*f(n)+1512*f(n+1)):

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

err := Vector(round(N)):
exy_lst := Vector(round(N)):
numerical_y1:=Vector(round(N)):

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

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

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


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

	par1:=x[0]=t[0],x[1]=t[1],x[2]=t[2]:
	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
		printf("%6.5f%45.30f%45.30f\n", 
		h*c,res[i],res[i+2]):
		
		numerical_y1[c] := res[i]:
		
		c:=c+1:
	end do:
	iny:=res[2]:
	ind:=res[4]:
	inx:=t[2]:
	for j from 0 to 2 do
		t[j]:=inx + j*h:
	end do:
end do:
v:=time() - st;
v/4;
printf("Maximum error is %.13g\n", max(err));
NFE=evalf((N/4*3)+1);
#get array of numerical and exact solutions for y1
numerical_array_y1 := [seq(numerical_y1[i], i = 1 .. N)]:
#exact_array_y1 := [seq(exy[i], i = 1 .. N)]:

#get array of time steps
time_t := [seq](step[i], i = 1 .. N):

#display graphs for y1
with(plots):
numerical_plot_y1 := plot(time_t, numerical_array_y1, style = [point], symbol = [asterisk],
				color = [blue,blue],symbolsize = 20, legend = ["TFIBF"]);

 Thank you, and best regards.

How can numbers be displayed inside the contour plot?

 restart;
with(plots);
contourplot(x*exp(-x^2 - y^2), x = -2 .. 2, y = -2 .. 2, axes = boxed);
like this

Dear Colleagues,

I wish to use plot3d to the attached code but always encoutered error. However, pointplot3d runs perfectly. Please I need your assistance in this regards.

Thank you all and best regards.K2_Problem_2_two_body_kepler_e=0.mw

Download Ex2_FDAEs.mw

I want to calculate the double integral of the following expression which includes sum of several Legendre polynomial terms, but the speed is so low. Any suggestion to speed up the calculation?

Download Soal.mw

Dear esteem Colleagues,

Please how do I modify the following two files (though similar) to get consistent errors? I am not sure where I made the mistake.

Any modifications would be appreciated.

Thank you all for your time and mentorship. Best regard

Biratu_Mapleprimes.mw

DDE_2_Mapleprime.mw

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
  q[r]:=(x^r-add(evalf(Int(x^r*q[R]*G,x=(-1)..1))*q[R]/evalf(Int(q[R]^2*G,x=(-1)..1)),R=0..r-1));
  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
      M(s+1,j+1):=add(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;
    end do;
    
    M(s+1,n+1):=add(coeff(a[n][s],x,k)*coeff(Koeffizienten,x,k),k=0..n);
    
    
  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;

Ergebnis:=add(eval(f,x=KnotenHnm[k])*Gewichte[k],k=1..nops([KnotenHnm]));

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

     a[s][j]:=c[j-1]*add(coeff(a[s][j-1],x,k-1)/c[k-1]*x^k,k=abs(s-j)+2..s+j)+add((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]*add(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]*add(coeff(a[s][j-2],x,k)*x^k,k=abs(s-j)+2..s+j-2);

     
    
  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?

Thank you in advance!

Hi,

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

Thanks 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

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?

Download ask_maple.mw

doubt_1.mw

Hi, I am trying to solve two simultaneous equations (for t1) they are as follows-

eq 1

i__m2(0) = (-(-b*`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`+`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`^2-2*c)*exp(0)*a*N^alpha/`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`^3+(-c*t__2^2*`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`^2+b*t__2*`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`^2+2*c*t__2*`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`-b*`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`+`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`^2-2*c)*exp(`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`*t__2)*a*N^alpha/`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`^3)*exp(0)

eq 2

i__m1(t__1) = ((-c*t__1^2*`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`^2+b*t__1*`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`^2+2*c*t__1*`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`-b*`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`+`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`^2-2*c)*exp(`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`*t__1)*a*N^alpha*(lambda-1)/`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`^3-(-b*`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`+`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`^2-2*c)*a*N^alpha*(lambda-1)/`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`^3)*exp(-`#msub(mi("&theta;",fontstyle = "normal"),mi("m"))`*t__1)

rhs(i__m2(0) = (-(-b*theta__m+theta__m^2-2*c)*exp(0)*a*N^alpha/theta__m^3+(-c*t__2^2*theta__m^2+b*t__2*theta__m^2+2*c*t__2*theta__m-b*theta__m+theta__m^2-2*c)*exp(theta__m*t__2)*a*N^alpha/theta__m^3)*exp(0)) = rhs(i__m1(t__1) = ((-c*t__1^2*theta__m^2+b*t__1*theta__m^2+2*c*t__1*theta__m-b*theta__m+theta__m^2-2*c)*exp(theta__m*t__1)*a*N^alpha*(lambda-1)/theta__m^3-(-b*theta__m+theta__m^2-2*c)*a*N^alpha*(lambda-1)/theta__m^3)*exp(-theta__m*t__1))

solve({-(-b*theta__m+theta__m^2-2*c)*a*N^alpha/theta__m^3+(-c*t__2^2*theta__m^2+b*t__2*theta__m^2+2*c*t__2*theta__m-b*theta__m+theta__m^2-2*c)*exp(theta__m*t__2)*a*N^alpha/theta__m^3 = ((-c*t__1^2*theta__m^2+b*t__1*theta__m^2+2*c*t__1*theta__m-b*theta__m+theta__m^2-2*c)*exp(theta__m*t__1)*a*N^alpha*(lambda-1)/theta__m^3-(-b*theta__m+theta__m^2-2*c)*a*N^alpha*(lambda-1)/theta__m^3)*exp(-theta__m*t__1)}, [t__1]);
Warning, solutions may have been lost
 

Can someone, please help. Thanks in advance.

Dear All.

I hope we are all staying safe.

Please I want to solve Sine Gordon Equation using a numerical method I constructed (non-classical), I need to compare the result of the method with the exact solution to generate the errors. However, since the exact solution has two variables, x, and t, I really don't know how to accommodate the two in my coding.

Can someone be of help in this regard?

Thank you and kind regards

Download Discretization_of_Sine_Gordon_Equation.mw

Download Sine_Gordon_Implementation_Trial.mw

When executing DEBUG within inline code (not within a procedure) the values displayed in successive debug windows (on clicking continue) are added to the end of my worksheet. How can the latter display be prevented?

Please I found out that the MatrixInverse on the assignment statement P3 does not run for about three days now. Please kindly help to simplify the code. Thank you and kind regards.

restart; omega := v/h;
r := a[0]+a[1]*x+a[2]*sinh(omega*x)+a[3]*cosh(omega*x)+a[4]*cos(omega*x)+a[5]*sin(omega*x);
b := diff(r, x);

c := eval(b, x = q) = f[n];
d := eval(r, x = q+3*h) = y[n+3]; e := eval(b, x = q+3*h) = f[n+3];
g := eval(b, x = q+2*h) = f[n+2];
i := eval(b, x = q+h) = f[n+1];
j := eval(b, x = q+4*h) = f[n+4];
k := solve({c, d, e, g, i, j}, {a[0], a[1], a[2], a[3], a[4], a[5]});
Warning,  computation interrupted
assign(k);
cf := r;
s4 := y[n+4] = simplify(eval(cf, x = q+4*h));
s3 := y[n+2] = simplify(eval(cf, x = q+2*h));
s2 := y[n+1] = simplify(eval(cf, x = q+h));
s1 := y[n] = simplify(eval(cf, x = q));

with(LinearAlgebra);
with(plots);
h := 1;
YN_1 := seq(y[n+k], k = 1 .. 4);
A1, a0 := GenerateMatrix([s1, s2, s3, s4], [YN_1]);
eval(A1);
YN := seq(y[n-k], k = 3 .. 0, -1);
A0, b1 := GenerateMatrix([s1, s2, s3, s4], [YN]);
eval(A0);
FN_1 := seq(f[n+k], k = 1 .. 4);
B1, b2 := GenerateMatrix([s1, s2, s3, s4], [FN_1]);
eval(B1);
FN := seq(f[n-k], k = 3 .. 0, -1);
B0, b3 := GenerateMatrix([s1, s2, s3, s4], [FN]);
eval(B0);
ScalarMultiply(R, A1)-A0;
det := Determinant(ScalarMultiply(R, A1)-A0);
P1 := A1-ScalarMultiply(B1, z);
P2 := combine(simplify(P1, size), trig);
P3 := MatrixInverse(P2);
P4 := A0-ScalarMultiply(B0, z);
P5 := MatrixMatrixMultiply(P3, P4);
P6 := Eigenvalues(P5);
f := P6[4];
T := unapply(f, z);
implicitplot(f, z = -5 .. 5, v = -5 .. 5, filled = true, grid = [5, 5], gridrefine = 8, labels = [z, v], coloring = [blue, white]);

Hello,

I want to sort the formulae to Psi and Beta, but I don't know how it works. I have tried it with sort, simplify, isolate, but that isn't what I'm searching.

It should looks like the simplier formula in the picture.

ab := (diff(Psii(t), t, t))*J-l[f]*(F[s, f, l]+F[s, f, r])+l[r]*(F[s, r, l]+F[s, r, r])-(1/2)*b[r]*(-F[s, r, l]*delta[l]+F[s, r, r]*delta[r]) = 0;
  / d  / d         \\                                   
  |--- |--- Psii(t)|| J - l[f] (F[s, f, l] + F[s, f, r])
  \ dt \ dt        //                                   

     + l[r] (F[s, r, l] + F[s, r, r])

       1                                                      
     - - b[r] (-F[s, r, l] delta[l] + F[s, r, r] delta[r]) = 0
       2                                                      
bc := (diff(betaa(t), t, t))*m*v*betaa(t)+F[s, r, l]*delta[l]+F[s, r, r]*delta[r]-(diff(Psii(t), t)) = 0;
    / d  / d          \\                                   
    |--- |--- betaa(t)|| m v betaa(t) + F[s, r, l] delta[l]
    \ dt \ dt         //                                   

                               / d         \    
       + F[s, r, r] delta[r] - |--- Psii(t)| = 0
                               \ dt        /    
cd := (diff(betaa(t), t))*m*v+F[s, r, l]+F[s, r, r]+F[s, f, l]+F[s, f, r]-(diff(Psii(t), t)) = 0;
   / d          \                                           
   |--- betaa(t)| m v + F[s, r, l] + F[s, r, r] + F[s, f, l]
   \ dt         /                                           

                     / d         \    
      + F[s, f, r] - |--- Psii(t)| = 0
                     \ dt        /    
F[s, f, l] := c[fl]*alpha[fl];
                        c[fl] alpha[fl]
F[s, f, r] := c[fr]*alpha[fr];
                        c[fr] alpha[fr]
F[s, r, l] := c[rl]*alpha[rl];
                        c[rl] alpha[rl]
F[s, r, r] := c[rr]*alpha[rr];
                        c[rr] alpha[rr]
alpha[fl] := (-v*betaa-l[f]*(diff(Psii(t), t)))/(-v+(1/2)*b[f]*(diff(Psii(t), t)));
                                 / d         \
                 -v betaa - l[f] |--- Psii(t)|
                                 \ dt        /
                 -----------------------------
                        1      / d         \  
                   -v + - b[f] |--- Psii(t)|  
                        2      \ dt        /  
alpha[fr] := (-v*betaa-l[f]*(diff(Psii(t), t)))/(v-(1/2)*b[f]*(diff(Psii(t), t)));
                                 / d         \
                 -v betaa - l[f] |--- Psii(t)|
                                 \ dt        /
                 -----------------------------
                       1      / d         \   
                   v - - b[f] |--- Psii(t)|   
                       2      \ dt        /   
alpha[rl] := delta[l]+(-v*betaa+l[r]*(diff(Psii(t), t)))/(-v+(1/2)*b[r]*(diff(Psii(t), t)));
                                       / d         \
                       -v betaa + l[r] |--- Psii(t)|
                                       \ dt        /
            delta[l] + -----------------------------
                              1      / d         \  
                         -v + - b[r] |--- Psii(t)|  
                              2      \ dt        /  
alpha[rr] := delta[r]+(-v*betaa+l[r]*(diff(Psii(t), t)))/(-v-(1/2)*b[r]*(diff(Psii(t), t)));
                                       / d         \
                       -v betaa + l[r] |--- Psii(t)|
                                       \ dt        /
            delta[r] + -----------------------------
                              1      / d         \  
                         -v - - b[r] |--- Psii(t)|  
                              2      \ dt        /  

ab;
                             /
                             |
/ d  / d         \\          |
|--- |--- Psii(t)|| J - l[f] |
\ dt \ dt        //          |
                             |
                             \

        /                / d         \\
  c[fl] |-v betaa - l[f] |--- Psii(t)||
        \                \ dt        //
  -------------------------------------
             1      / d         \      
        -v + - b[f] |--- Psii(t)|      
             2      \ dt        /      

           /                / d         \\\        /      /      
     c[fr] |-v betaa - l[f] |--- Psii(t)|||        |      |      
           \                \ dt        //|        |      |      
   + -------------------------------------| + l[r] |c[rl] |delta[
               1      / d         \       |        |      |      
           v - - b[f] |--- Psii(t)|       |        |      |      
               2      \ dt        /       /        \      \      

                       / d         \\
       -v betaa + l[r] |--- Psii(t)||
                       \ dt        /|
  l] + -----------------------------|
              1      / d         \  |
         -v + - b[r] |--- Psii(t)|  |
              2      \ dt        /  /

           /                           / d         \\\          /
           |           -v betaa + l[r] |--- Psii(t)|||          |
           |                           \ dt        /||   1      |
   + c[rr] |delta[r] + -----------------------------|| - - b[r] |
           |                  1      / d         \  ||   2      |
           |             -v - - b[r] |--- Psii(t)|  ||          |
           \                  2      \ dt        /  //          \
       /                           / d         \\         
       |           -v betaa + l[r] |--- Psii(t)||         
       |                           \ dt        /|         
-c[rl] |delta[l] + -----------------------------| delta[l]
       |                  1      / d         \  |         
       |             -v + - b[r] |--- Psii(t)|  |         
       \                  2      \ dt        /  /         

           /                           / d         \\         \   
           |           -v betaa + l[r] |--- Psii(t)||         |   
           |                           \ dt        /|         |   
   + c[rr] |delta[r] + -----------------------------| delta[r]| = 
           |                  1      / d         \  |         |   
           |             -v - - b[r] |--- Psii(t)|  |         |   
           \                  2      \ dt        /  /         /   

  0
bc;
 / d  / d          \\             
 |--- |--- betaa(t)|| m v betaa(t)
 \ dt \ dt         //             

            /                           / d         \\         
            |           -v betaa + l[r] |--- Psii(t)||         
            |                           \ dt        /|         
    + c[rl] |delta[l] + -----------------------------| delta[l]
            |                  1      / d         \  |         
            |             -v + - b[r] |--- Psii(t)|  |         
            \                  2      \ dt        /  /         

            /                           / d         \\         
            |           -v betaa + l[r] |--- Psii(t)||         
            |                           \ dt        /|         
    + c[rr] |delta[r] + -----------------------------| delta[r]
            |                  1      / d         \  |         
            |             -v - - b[r] |--- Psii(t)|  |         
            \                  2      \ dt        /  /         

      / d         \    
    - |--- Psii(t)| = 0
      \ dt        /    
cd;
 / d          \    
 |--- betaa(t)| m v
 \ dt         /    

            /                           / d         \\
            |           -v betaa + l[r] |--- Psii(t)||
            |                           \ dt        /|
    + c[rl] |delta[l] + -----------------------------|
            |                  1      / d         \  |
            |             -v + - b[r] |--- Psii(t)|  |
            \                  2      \ dt        /  /

            /                           / d         \\
            |           -v betaa + l[r] |--- Psii(t)||
            |                           \ dt        /|
    + c[rr] |delta[r] + -----------------------------|
            |                  1      / d         \  |
            |             -v - - b[r] |--- Psii(t)|  |
            \                  2      \ dt        /  /

            /                / d         \\
      c[fl] |-v betaa - l[f] |--- Psii(t)||
            \                \ dt        //
    + -------------------------------------
                 1      / d         \      
            -v + - b[f] |--- Psii(t)|      
                 2      \ dt        /      

            /                / d         \\                    
      c[fr] |-v betaa - l[f] |--- Psii(t)||                    
            \                \ dt        //   / d         \    
    + ------------------------------------- - |--- Psii(t)| = 0
                1      / d         \          \ dt        /    
            v - - b[f] |--- Psii(t)|                           
                2      \ dt        /