Maple Codes

> Taylor3D:=proc (f::procedure, a::realcons, b::realcons,n::posint, rx::equation,ry::equation,rz::equation)
local t, i,pltf,disp,r,p;
readlib(mtaylor);
pltf:= plot3d(f(x,y),rx,ry,style=patch);
for i to n do
r(i):=mtaylor(f(x,y),[x = a,y=b],i+1);
t(i):=plots[textplot3d]([2,2,0.5,convert(Degree=i,string)],font=[TIMES,BOLD,14],color=red):
p(i) := plots[display](t(i),plot3d(r(i),rx,ry,style=patch,color =
COLOR(RGB,1/1000000000000*rand(),1/1000000000000*rand(),1/1000000000000*rand())));
disp(i):=plots[display](p(i),pltf);
od;
plots[display]([seq(disp(i),i = 1 .. n)],view = [rhs(rx), rhs(ry),rhs(rz)],axes=frame,insequence = true,args[8..nargs]);
end:

> Pendulum:=proc(L,theta0,c)
local disk1,disk2,disk3,plt,p,background1,background2,deq,init,alpha,line,sol,plts,counter,i,pltdisk,
period,counterback,hscale,theta01,tk;
with(plots):with(plottools):
deq:=diff(theta(t),t$2)+c*diff(theta(t),t)+9.8/L*sin(theta(t))=0;
print(deq);
theta01:=theta0/180*Pi;
init:=theta(0)=theta01, D(theta)(0)=0;
sol:=dsolve({deq,init},theta(t),numeric,startinit=true);
alpha:=u->subs(sol(u),theta(t));
period:=4*fsolve('alpha(t)',t=0..2);
tk:=2;
hscale:=tk*period;
if c=0 then print(Period) else print(Pseudoperiod) fi;
print(`T`=period);
counter:=t->textplot([L/2,0.15*L,convert([evalf(t,2),evalf(alpha(t)/Pi*180,4)],string)],
color=brown,font=[TIMES,BOLD,14]);
disk1:=disk([0,0],L/25,color=blue);
disk2:=t->translate(disk1,t,alpha(t));
plts:=plot(['alpha(t)'+L,L],t=0..tk*period,color=[aquamarine,orange],thickness=[2,1]);
pltdisk:=t->translate(disk2(t),0,L);
background1:=pieslice([0,0],21/20*L,-(Pi/2+theta01)..-(Pi/2-theta01),color=grey):
background2:=rectangle([-21/20*L,-21.5/20*L],[21/20*L,0.26*L],color=yellow):
counterback:=rectangle([0.1,0.01*L],[1.03*L,0.24*L],color=cyan):
line:=plot([[0,0],[0,-L]],thickness=2,color=red);
disk3:=disk([0,-L],L/20,color=blue);
plt:=display(line,disk3);
p:=t->display(rotate(plt,-alpha(t)),plts,pltdisk(t),counterback,background1,background2,counter(t),
view=[-1.1*L..1.1*hscale,-1.1*L..L+alpha(0)+L/25]);
display(seq(p(i/10),i=0..10*hscale),insequence=true,args[4..nargs]);
end:

> SpringMassCouplet:=proc(L::list,init::list,tr)
local spring1,spring2,spring3,spring4,mass1,mass2,mass3,deq1,deq2,deq3,init2,
sol,xk1,xk2,xk3,plt1,plt2,plt, pltxk,base,M,K,A, a, yt,xs,AA,
disk1,disk2,disk3,initpos2,initpos3,spring,m1,m2,m3,k1,k2,k3,k4,s3,s4,eq;
with(plots):with(plottools):with(linalg):
spring:=proc(x1,x2,n,c)
local p1,p2,p3,p4,pn_1,pn,p;
p1:=[x1,0];
p2:=[x1+0.25,0];
p3:=[x2-0.25,0];
p4:=[x2,0];
p:=i->[x1+0.25+(x2-x1-0.5)/n*i,(-1)^(i+1)*0.5];
plot([p1,p2,seq(p(i),i=1..n-1),p3,p4],thickness=2,color=c);
end:
a:=2; yt:=6; xs:=tr;
m1:=L[1,1]: m2:=L[2,1]:
k1:=L[1,2]: k2:=L[2,2]:
base:=plot([[0,y,y=-1.1..1.6],[21.9,y,y=-1.1..1.6],[t,-1.1,t=0..22]],color=brown,thickness=4);
spring1:=x1->spring(0,x1,12,aquamarine);
spring2:=(x1,x2)->spring(x1+2,x2-2,12,aquamarine);
spring3:=(x2,x3)->spring(x2,x3-2,12,aquamarine);
spring4:=x3->spring(x3,22,12,aquamarine);
mass1:=x1->rectangle([x1,-1],[x1+2,1.6],color=red):
mass2:=x2->rectangle([x2-2,-1],[x2,1.6],color=blue):
mass3:=x3->rectangle([x3-2,-1],[x3,1.6],color=grey):
if nops(L) =4 then
m3:=L[3,1]; k3:=L[3,2]; k4:=L[4];
s4:=i->spring4(xk3(xs/40*i)+16):
fi;
if nops(L)=3 and nops(L[3])=2 then
m3:=L[3,1]; k3:=L[3,2]; k4:=0;
fi;
if nops(L)>2 and nops(L[3])=2 then
deq1:=m1*diff(x1(t),t$2)+k1*x1(t)-k2*(x2(t)-x1(t))=0;
deq2:=m2*diff(x2(t),t$2)+k2*(x2(t)-x1(t))-k3*(x3(t)-x2(t))=0;
deq3:=m3*diff(x3(t),t$2)+k3*(x3(t)-x2(t))+k4*x3(t)=0;
print(deq1,deq2,deq3);
init2:=x1(0)=init[1], x2(0)=init[2], x3(0)=init[3], D(x1)(0)=init[4],
D(x2)(0)=init[5],D(x3)(0)=init[6];
sol:=dsolve({deq1,deq2,deq3,init2},{x1(t),x2(t),x3(t)},method=laplace);
sol:=combine(simplify(sol));
xk1:=unapply(subs(sol,x1(t)),t);
xk2:=unapply(subs(sol,x2(t)),t);
xk3:=unapply(subs(sol,x3(t)),t);
pltxk:=plot([xk1(t)+2,xk2(t)+a,xk3(t)+a,a],t=0..22,color=[red,blue,grey,black],
numpoints=400);
disk1:=t->translate(disk([0,0],0.2,color=red),t,xk1(t)+a);
disk2:=t->translate(disk([0,0],0.2,color=blue),t,xk2(t)+a);
disk3:=t->translate(disk([0,0],0.2,color=grey),t,xk3(t)+a);
initpos3:=plot([[5,y,y=-1..2],[10,y,y=-1..2],[15,y,y=-1..2]],color=[red,blue,grey],
thickness=2):
if nops(L)=3 and nops(L[3])=2 then
plt1:=i->translate(display(spring1(xk1(xs/40*i)+4),mass1(xk1(xs/40*i)+4),
spring2(xk1(xs/40*i)+4,xk2(xs/40*i)+11),
mass2(xk2(xs/40*i)+11),spring3(xk2(xs/40*i)+11,xk3(xs/40*i)+16),
mass3(xk3(xs/40*i)+16),base,initpos3),0,-5);
fi;
if nops(L)=4 and nops(L[4])=1 then
plt1:=i->translate(display(spring1(xk1(xs/40*i)+4),mass1(xk1(xs/40*i)+4),
spring2(xk1(xs/40*i)+4,xk2(xs/40*i)+11),
mass2(xk2(xs/40*i)+11),spring3(xk2(xs/40*i)+11,xk3(xs/40*i)+16),
mass3(xk3(xs/40*i)+16),s4(i),base,initpos3),0,-5);
fi;
M:=(m1,m2,m3)->diag(m1,m2,m3);
#print(`mass matrix`);
#print(M(m[1],m[2],m[3])=M(m1,m2,m3));
if nops(L)=4 then
K:=(k1,k2,k3,k4)->matrix(3,3,[-(k1+k2),k2,0,k2,-(k2+k3),k3,0,k3,-(k3+k4)]);
A:=(m1,m2,m3,k1,k2,k3,k4)->evalm(inverse(M(m1,m2,m3))&*K(k1,k2,k3,k4));
#print(`stiffness matrix`);
#print(K(k[1],k[2],k[3],k[4])=K(k1,k2,k3,k4));
print(`coefficient matrix`);
print(inverse(M(m[1],m[2],m[3]))*K(k[1],k[2],k[3],k[4])=A(m1,m2,m3,k1,k2,k3,k4));
print(eigenvalues);
AA:=A(m1,m2,m3,k1,k2,k3,k4);
fi;
if nops(L)=3 and nops(L[3])=2 then
K:=(k1,k2,k3)->matrix(3,3,[-(k1+k2),k2,0,k2,-(k2+k3),k3,0,k3,-k3]);
A:=(m1,m2,m3,k1,k2,k3)->evalm(inverse(M(m1,m2,m3))&*K(k1,k2,k3));
#print(`stiffness matrix`);
#print(K(k[1],k[2],k[3])=K(k1,k2,k3));
print(`coefficient matrix`);
print(inverse(M(m[1],m[2],m[3]))*K(k[1],k[2],k[3])=A(m1,m2,m3,k1,k2,k3));
print(eigenvalues);
AA:=A(m1,m2,m3,k1,k2,k3);
fi;
eq:= [-omega[1]^2,-omega[2]^2,-omega[3]^2]=
fnormal(evalf([eigenvals(AA)]));
print(eq);
print(`natural frequencies`);
print(map(sqrt,-lhs(eq),symbolic)=map(sqrt,-rhs(eq)));
print(x1(t)=xk1(t),x2(t)=xk2(t),x3(t)=xk3(t));
plt:=i->display(plt1(i),pltxk,disk1(xs/40*i),disk2(xs/40*i),disk3(xs/40*i),view=[0..22,-6.2..yt]);
print(display(seq(plt(i),i=0..40),insequence=true,args[4..nargs]));
fi;
if nops(L)>=3 and nops(L[3])=2 then RETURN();
else;
if nops(L)=2 and nops(L[2])=2 then
k3:=0;
fi;
if nops(L)=3 then k3:=L[3]; fi;
init2:=x1(0)=init[1], x2(0)=init[2], D(x1)(0)=init[3],D(x2)(0)=init[4];
deq1:=m1*diff(x1(t),t$2)+k1*x1(t)-k2*(x2(t)-x1(t))=0;
deq2:=m2*diff(x2(t),t$2)+k2*(x2(t)-x1(t))+k3*x2(t)=0;
print(deq1,deq2);
sol:=dsolve({deq1,deq2,init2},{x1(t),x2(t)},method=laplace);
sol:=combine(simplify(sol));
xk1:=unapply(subs(sol,x1(t)),t);
xk2:=unapply(subs(sol,x2(t)),t);
pltxk:=plot([xk1(t)+a,xk2(t)+a,a],t=0..22,color=[red,blue,black],numpoints=400);
disk1:=t->translate(disk([0,0],0.2,color=red),t,xk1(t)+a);
disk2:=t->translate(disk([0,0],0.2,color=blue),t,xk2(t)+a);
initpos2:=plot([[7,y,y=-1..2],[14,y,y=-1..2]],color=[red,blue],thickness=2):
spring3:=x2->spring(x2,22,12,aquamarine);
if nops(L)=2 and nops(L[2])=2 then
plt1:=i->translate(display(spring1(xk1(xs/40*i)+6),mass1(xk1(xs/40*i)+6),
spring2(xk1(xs/40*i)+6,xk2(xs/40*i)+15),
mass2(xk2(xs/40*i)+15),base,initpos2),0,-5);
fi;
if nops(L)=3 and nops(L[3])=1 then
k3:=L[3]; s3:=i->spring3(xk2(xs/40*i)+15);
plt1:=i->translate(display(spring1(xk1(xs/40*i)+6),mass1(xk1(xs/40*i)+6),
spring2(xk1(xs/40*i)+6,xk2(xs/40*i)+15),
mass2(xk2(xs/40*i)+15),s3(i),base,initpos2),0,-5);
fi;
if nops(L)=3 and nops(L[3])=1 then
M:=(m1,m2)->diag(m1,m2);
#print(`mass matrix`);
#print(M(m[1],m[2])=M(m1,m2));
K:=(k1,k2,k3)->matrix(2,2,[-(k1+k2),k2,k2,-(k2+k3)]);
#print(`stiffness matrix`);
#print(K(k[1],k[2],k[3])=K(k1,k2,k3));
A:=(m1,m2,k1,k2,k3)->evalm(inverse(M(m1,m2))&*K(k1,k2,k3));
print(`coefficient matrix`);
print(inverse(M(m[1],m[2]))*K(k[1],k[2],k[3])=A(m1,m2,k1,k2,k3));
print(eigenvalues);
eq:= [-omega[1]^2,-omega[2]^2]=fnormal(evalf([eigenvals(A(m1,m2,k1,k2,k3))]));
print(eq);
print(`natural frequencies`);
print(map(sqrt,-lhs(eq),symbolic)=map(sqrt,-rhs(eq)));
print(x1(t)=xk1(t),x2(t)=xk2(t));
plt:=i->display(plt1(i),pltxk,disk1(xs/40*i),disk2(xs/40*i),view=[0..22,-6.2..yt]);
print(display(seq(plt(i),i=0..40),insequence=true,args[4..nargs]));
fi;
if nops(L)=2 and nops(L[2])=2 then
M:=(m1,m2)->diag(m1,m2);
#print(`mass matrix`);
#print(M(m[1],m[2])=M(m1,m2));
K:=(k1,k2)->matrix(2,2,[-(k1+k2),k2,k2,-k2]);
#print(`stiffness matrix`);
#print(K(k[1],k[2])=K(k1,k2));
A:=(m1,m2,k1,k2)->evalm(inverse(M(m1,m2))&*K(k1,k2));
print(`coefficient matrix`);
print(inverse(M(m[1],m[2]))*K(k[1],k[2])=A(m1,m2,k1,k2));
print(eigenvalues);
eq:= [-omega[1]^2,-omega[2]^2]=fnormal(evalf([eigenvals(A(m1,m2,k1,k2))]));
print(eq);
print(`natural frequencies`);
print(map(sqrt,-lhs(eq),symbolic)=map(sqrt,-rhs(eq)));
print(x1(t)=xk1(t),x2(t)=xk2(t));
plt:=i->display(plt1(i),pltxk,disk1(xs/40*i),disk2(xs/40*i),view=[0..22,-6.2..yt]);
display(seq(plt(i),i=0..40),insequence=true,args[4..nargs]);
fi;
fi;
end:

> WavePulse:=proc(m1,m2 )
local finn,ft,pulsi,pulsr,pulst,ft0,fis,plt,r,tau,fr,tykv,tykh,background;
with(plots):with(plottools):
background:=rectangle([-10,-4],[10,4],color=yellow):
r:=(sqrt(m1)-sqrt(m2))/(sqrt(m1)+sqrt(m2));
tau:=2*sqrt(m1)/(sqrt(m1)+sqrt(m2));
finn:=x->piecewise(x<-2.4,0.02,x<2.4,2.5*exp(-x^2),0.02):
fr:=x->piecewise(x<-2.4,0.02,x<2.4,r*2.5*exp(-x^2),0.02):
ft:=x->piecewise(x<0,0,x<4.8,tau*2.5*exp(-(x-2.4)^2),0.02):
if m1<m2 then tykv:=2; tykh:=4;
elif m1=m2 then tykv:=2; tykh:=2;
else tykv:=4;tykh:=2;
fi;
ft0:=plot(0,x=0..10,color=blue,thickness=tykh);
fis:=plot(0,x=-10..0,color=red,thickness=tykv);
pulsi:=t->plot(finn(x-t+2.4),x=-10..0,color=red,scaling=constrained,thickness=tykv):
pulsr:=t->plot(fr(x+t),x=-10..0,color=red,scaling=constrained,thickness=tykv):
pulst:=t->plot(ft(x-t),x=0..10,color=blue,scaling=constrained,thickness=tykh):
plt:=t->if t<0 then display(pulsi(t),ft0,background); else display(pulst(t),pulsr(t+2.4),background); fi;
display(seq(plt(t),t=-10..10),insequence=true,scaling=constrained,args[3..nargs]);
end:

> Convolution:=proc(f::procedure,g::procedure)
local r1,pltf,pltg,plt2,h,plth,plt,text,disk1,plth2;
with(plots):with(plottools):
alias(u=Heaviside):
r1:=display(rectangle([0,0],[1,1],color=red)):
pltf:=plot(f(t),t=-4..2,color=red,filled=true):
pltg:=plot(g(t)*u(t),t=-7..7,y=-0.2..1.2,color=green,filled=true):
plt2:=tk->if tk<-1 then display(translate(r1,tk,0),pltg) else display(pltg,translate(r1,tk,0)) fi:
h:=t->Int('f(tau)*g(t-tau)',tau=0..t);
print(h(t)=value(h(t)));
plth:=plot(value(h(t)),t=0..6,color=blue,numpoints=400,title=`Convolution product`);
text:=tk->textplot([0.5,1.1,convert(Overlapping_Area=evalf(h(tk+1)),string)],
font=[TIMES,BOLD,14],color=brown);
disk1:=tk->translate(ellipse([0,0],0.2,0.04,color=green,filled=true),tk+1,value(h(tk+1)));
plth2:=tk->display(plth,disk1(tk));
plt:=tk->display(plt2(tk),text(tk),translate(plth2(tk),-6,0));
display(seq(plt(-2+7/40*i),i=0..40),insequence=true,tickmarks=[4,0]);
end:

> PointIteration:=proc(A::matrix,LX0::list,N::integer,Style)
local n,pltL1,pltL2,pltL3,plt,ev,i,L,X,points,xmax,xmin,ymax,ymin,background,hscale,vscale;
with(linalg);with(plottools):
ev:=eigenvals(A);
for i from 1 to nops(LX0) do
X(i,0):=LX0[i];
for n from 1 to N do
X(i,n):=convert(evalm(A&*X(i,n-1)),list);
od:
od;
points:=seq(seq(X(i,n),i=1..nops(LX0)),n=0..N);
xmin:=min(op(map2(op,1,[points])));
xmax:=max(op(map2(op,1,[points])));
ymin:=min(op(map2(op,2,[points])));
ymax:=max(op(map2(op,2,[points])));
hscale:=xmax-xmin;
vscale:=ymax-ymin:
background:=rectangle([xmin-0.1*hscale,ymin-0.1*vscale],[xmax+0.1*hscale,ymax+0.1*vscale],
color=yellow);
L:=(k,n)->[seq(X(k,i),i=0..n)];
pltL1:=(k,n)->plot([L(k,n)],color=blue,style=Style,symbol=circle);
pltL2:=(k,n)->plot([L(k,n)],color=blue,style=point,symbol=circle);
pltL3:=(k,n)->plot([X(k,n)],color=red,style=point,symbol=circle);
plt:=n->display(seq(plots[display](pltL3(k,n),pltL1(k,n),pltL2(k,n),background),k=1..nops(LX0)));
print(A=evalm(A));
print(`Eigenvalues`);
print([lambda[1],lambda[2]]=[ev]);
plots[display](seq(plt(n),n=0..N),insequence=true,args[5..nargs]);
end:

>

>