@helix
Here is the code.
restart;
Code Packages
with(plots):
with(CodeGeneration):
with(LinearAlgebra):
Important Numbers and Functions
1 Square Numbers and Functions
phi:=proc(n,m,l1,l2)
RETURN((2/sqrt(l1*l2))*sin(((n*Pi)/l1)*x)*sin(((m*Pi)/l2)*y));
end proc:
lamda:=proc(n,m,l1,l2)
RETURN(((n*Pi)/l1)^2 + ((m*Pi)/l2)^2);
end proc:
omega:=proc(n,m,l1,l2)
RETURN(sqrt(lamda(n,m,l1,l2)));
end proc:
wave:=proc(n,m,l1,l2)
RETURN(phi(n,m,l1,l2) * cos(omega(n,m,l1,l2) * t));
end proc:
wave2:=proc(n,m,l1,l2)
RETURN(phi(n,m,l1,l2) * sin(omega(n,m,l1,l2) * t));
end proc:
2 Squares Numbers and Functions
HeavyStepLeft:= piecewise(x<0,1);
HeavyStepRight:= piecewise(x>=0,1);
aPrime:= proc(n1,n2,a,m,l2)
local AP;
AP:= sqrt(((n2^2/n1^2)*a^2) + ((m^2/l2^2)*Pi^2)*((n2^2/n1^2) -1 ));
RETURN(AP);
end proc:
interval_array_ten:= Array([[1.25,1.75,2,2.25,2.75,3,3.5,3.75,4,4.5],[1.25,1.75,2,2.5,2.75,3.25,3.5,4,4.25,4.75],[1.25,2,2.5,2.75,3.25,3.75,4,4.25,4.75,5],[1.75,2.25,2.75,3.25,3.5,4,4.25,4.75,5,5.5],[1.25,2,2.5,3,3.5,4,4.25,4.75,5,5.5],[1.25,2.25,2.75,3.25,3.75,4.25,4.5,5,5.5,6],[1.75,2.5,3,3.5,4,4.5,5,5.25,5.75,6.25],[1.25,2,2.5,3.25,3.75,4.25,4.75,5.25,5.75,6.25],[1.25,2.25,2.75,3.5,4,4.5,5,5.5,6,6.5],[1.5,2.5,3,3.75,4.25,4.75,5.25,5.75,6.25,6.75]]);
A_norm:= proc(n1,n2,l1,l2,a,m)
local A;
if(sin(aPrime(n1,n2,a,m,l2)*l1) = 0)
then
A:= aPrime(n1,n2,a,m,l2)*cos(aPrime(n1,n2,a,m,l2)* l1) * (1/(sqrt(((aPrime(n1,n2,a,m,l2))^2 * cos(aPrime(n1,n2,a,m,l2) * l1)^2)* (l2/2) * (n1^2) * ((-sin(2*a*l1) + 2*a*l1)/(4*a)) + a^2 *(cos(a*l1))^2 *(l2/2) * (n2^2) * ((-sin(2*aPrime(n1,n2,a,m,l2)*l1) + 2*aPrime(n1,n2,a,m,l2)*l1)/(4*aPrime(n1,n2,a,m,l2))))));
else
A:= sin(aPrime(n1,n2,a,m,l2)* l1) * (1/(sqrt( (sin(aPrime(n1,n2,a,m,l2)*l1))^2* (l2/2) * (n1^2) * ((-sin(2*a*l1) + 2*a*l1)/(4*a)) + (sin(a*l1))^2*(l2/2) * (n2^2) * ((-sin(2*aPrime(n1,n2,a,m,l2)*l1) + 2*aPrime(n1,n2,a,m,l2)*l1)/(4*aPrime(n1,n2,a,m,l2))))));
fi;
RETURN(evalf(A));
end proc:
APrime_norm:= proc(n1,n2,l1,l2,a,m)
local A;
if(sin(a*l1) = 0)
then
A:= -a*cos(a*l1) * (1/(sqrt(((aPrime(n1,n2,a,m,l2))^2 * cos(aPrime(n1,n2,a,m,l2) * l1)^2)* (l2/2) * (n1^2) * ((-sin(2*a*l1) + 2*a*l1)/(4*a)) + a^2 *(cos(a*l1)^2)*(l2/2) * (n2^2) * ((-sin(2*aPrime(n1,n2,a,m,l2)*l1) + 2*aPrime(n1,n2,a,m,l2)*l1)/(4*aPrime(n1,n2,a,m,l2))))));
else
A:= sin(a*l1) * (1/(sqrt((sin(aPrime(n1,n2,a,m,l2) * l1)^2)* (l2/2) * (n1^2) * ((-sin(2*a*l1) + 2*a*l1)/(4*a)) + (sin(a*l1))^2 *(l2/2) * (n2^2) * ((-sin(2*aPrime(n1,n2,a,m,l2)*l1) + 2*aPrime(n1,n2,a,m,l2)*l1)/(4*aPrime(n1,n2,a,m,l2))))));
fi;
RETURN(evalf(A));
end proc:
lambda_Two_squares:= proc(n_value,a_value,m,l2)
RETURN((1/n_value^2)*(a_value^2 + (m^2/l2^2)*Pi^2));
end proc:
omega_Two_squares:= proc(n_value,a_value,m,l2)
RETURN(sqrt(lambda_Two_squares(n_value,a_value,m,l2)));
end proc:
Phi_left:= proc(n1,n2,l1,l2,a,m)
local A;
A:= A_norm(n1,n2,l1,l2,a,m);
RETURN(A*sin(a*(x+l1))*sin(((m*Pi)/l2)*y));
end proc:
Phi_right:= proc(n1,n2,l1,l2,a,m)
local A;
local ap;
A:= APrime_norm(n1,n2,l1,l2,a,m);
ap:= aPrime(n1,n2,a,m,l2);
RETURN((A*sin(ap*(l1-x))*sin(((m*Pi)/l2)*y)));
end proc:
PhiFunc:= proc(n1,n2,l1,l2,a,m)
RETURN(Phi_left(n1,n2,l1,l2,a,m) *HeavyStepLeft + Phi_right(n1,n2,l1,l2,a,m) *HeavyStepRight);
end proc:
n_of_x:= piecewise(x<0,n1,x>0,n2);
wave_Two_squares:= proc(n1,n2,l1,l2,a,m)
RETURN(PhiFunc(n1,n2,l1,l2,a,m)* cos(omega_Two_squares(n1,a,m,l2) * t));
end proc:
wave_Two_squares2:= proc(n1,n2,l1,l2,a,m)
RETURN(PhiFunc(n1,n2,l1,l2,a,m)* sin(omega_Two_squares(n1,a,m,l2) * t));
end proc:
Both
FMaker := proc(x,y,Vx,Vy,Sigma)
local f;
f:= exp(-Sigma*((x-(Pi/2))^2 + ((y-(Pi/2))^2)))*cos(Vx*x + Vy*y);
RETURN (f);
end proc:
GMaker := proc(f, Vx, Vy, Sigma)
local g;
local DX;
local DY;
local ExPart := (simplify(f/exp(-Sigma*(y-(Pi/2))^2)))/cos(Vx*x + Vy*y);
local EyPart := (simplify(f/exp(-Sigma*(x-(Pi/2))^2)))/cos(Vx*x + Vy*y);
DX := ((EyPart*cos(Vy*y))*diff(ExPart*cos(Vx*x),x)-((EyPart*sin(Vy*y))*diff(ExPart*sin(Vx*x),x)));
DY := ((ExPart*cos(Vx*x))*diff(EyPart*cos(Vy*y),y)-((ExPart*sin(Vx*x))*diff(EyPart*sin(Vy*y),y)));
g:= Vx*DX + Vy*DY;
RETURN (g);
end proc:
1 Square Functions
CoeffFinder:=proc(f,Sigma,sum,l1,l2,Vx,Vy)
local n;
local m;
local F;
local SinPartx;
local SinParty;
local C := 2/sqrt(l1*l2);
local ExPart := (simplify(f/exp(-Sigma*(y-(Pi/2))^2)))/cos(Vx*x + Vy*y);
local EyPart := (simplify(f/exp(-Sigma*(x-(Pi/2))^2)))/cos(Vx*x + Vy*y);
F:=Array(1..sum,1..sum);
for n from 1 by 1 to sum
do
for m from 1 by 1 to sum
do
SinPartx:= sin(((n*Pi)/l1)*x);
SinParty:= sin(((m*Pi)/l2)*y);
F[n,m]:= evalf(C*((evalf(int(ExPart*SinPartx*cos(Vx*x),x=0..l1,numeric=true))*evalf((int(EyPart*SinParty*cos(Vy*y),y=0..l2,numeric=true)))-(evalf((int(ExPart*SinPartx*sin(Vx*x),x=0..l1,numeric=true)))*evalf((int(EyPart*SinParty*sin(Vy*y),y=0..l2,numeric=true)))))));
end do;
end do;
RETURN(F);
end proc:
CoeffFinder2:=proc(g,sum,l1,l2,Vx,Vy)
local n;
local m;
local G;
G:=Array(1..sum,1..sum);
for n from 1 by 1 to sum
do
for m from 1 by 1 to sum
do
G[n,m]:= (evalf(int(int(g*phi(n,m,l1,l2),x=0..l1,numeric=true),y=0..l2,numeric=true))/omega(n,m,l1,l2));
end do;
end do;
RETURN(G);
end proc:
WaveMaker:=proc(F,G,sum,l1,l2)
local n;
local m;
local totalWave;
totalWave := 0;
for n from 1 by 1 to sum
do
for m from 1 by 1 to sum
do
totalWave := totalWave + (F[n,m] * wave(n,m,l1,l2)) + (G[n,m] * wave2(n,m,l1,l2));
end do;
end do;
RETURN(totalWave);
end proc:
AnimationMaker:= proc(Sigma,Vx,Vy,Nodes,l1,l2,Dimensional)
local f;
local g;
local F;
local G;
f:= FMaker(x,y,Vx,Vy,Sigma);
g:= GMaker(f,Vx,Vy,Sigma);
F:=CoeffFinder(f,Sigma,Nodes,l1,l2,Vx,Vy);
G:=CoeffFinder2(g,Nodes,l1,l2,Vx,Vy);
if(Dimensional)
then
(RETURN(animate(plot3d,[WaveMaker(F,G,Nodes,l1,l2),x=0..l1,y=0..l2],t=0..10,axes=frame,frames=150, shading = zhue,scaling=constrained)));
else
(RETURN(animate(densityplot,[WaveMaker(F,G,Nodes,l1,l2),x=0..l1,y=0..l2,contrast=0.5],t=0..10,frames=150,scaling=constrained)));
end if;
end proc:
1 Square Output
#AnimationMaker(25,1,2,10,Pi,Pi,false);
#Check initial function against f.
(*
Sigma:=10;
Vx:=1;
Vy:=2;
l1:=Pi;
l2:=Pi;
Nodes:=10;
f:= FMaker(x,y,Vx,Vy,Sigma);
g:= GMaker(f,Vx,Vy,Sigma);
F:=CoeffFinder(f,Sigma,Nodes,l1,l2,Vx,Vy);
G:=CoeffFinder2(g,Nodes,l1,l2,Vx,Vy);
Waveee:=WaveMaker(F,G,Nodes,l1,l2):
plot3d(f,x=0..l1,y=0..l2);
t:=0;
plot3d(Waveee,x=0..l1,y=0..l2);
plot3d((Waveee-f),x=0..l1,y=0..l2,scaling=constrained);
*)
2 Squares with Boundary Functions
a_root_finder:= proc(n1,n2,m,l1,l2,interval1,interval2)
local M;
local a;
M:= Matrix([[sin(a*l1),-sin(aPrime(n1,n2,a,m,l2)*l1)],[a*cos(a*l1),aPrime(n1,n2,a,m,l2)*cos(aPrime(n1,n2,a,m,l2)*l1)]]);
RETURN(fsolve(Determinant(M) = 0, a=interval1..interval2));
end proc:
a_array_assembler:= proc(n1,n2,l1,l2,no_of_roots,no_of_ms,interval_array)
local i;
local m_local;
local start_point;
local a_Array;
local interval1;
local interval2;
start_point:= 0.1;
a_Array:= Array(1..no_of_roots,1..no_of_roots);
for m_local from 1 by 1 to no_of_ms
do
start_point := 0.1;
for i from 1 by 1 to no_of_roots
do
interval1:= start_point;
interval2:= interval_array[m_local,i];
a_Array[m_local,i] := a_root_finder(n1,n2,m_local,l1,l2,interval1,interval2);
start_point := interval2;
od;
od;
RETURN(a_Array);
end proc:
#This generates the values of 'a' which will be used in the left and right Eigenfunctions. Each row corresponds to m=1..10.
a_array:= a_array_assembler(1,2,Pi,Pi,10,10,interval_array_ten);
Two_squares_CoeffFinder:= proc(f,nodes,n1,n2,l1,l2, a_Array)
local F;
local count_a;
local m;
local a;
F:= Array(1..nodes,1..nodes);
for m from 1 by 1 to nodes
do
for count_a from 1 by 1 to nodes
do
a:= a_array[m,count_a];
#F[m,count_a]:= ((n1^2)*(int(int((Phi_left(n1,n2,l1,l2,a,m))*f,x=-l1..0,numeric=true),y=0..l2,numeric=true))) + ((n2^2)*(int(int((Phi_right(n1,n2,l1,l2,a,m))*f,x=0..l1,numeric=true),y=0..l2,numeric=true)));
F[m,count_a]:= int(int((n_of_x)^2 * PhiFunc(n1,n2,l1,l2,a,m) * f,x=-l1..l1,numeric=true),y=0..l1,numeric=true);
od;
od;
RETURN(F);
end proc:
Two_squares_CoeffFinder2:= proc(g,nodes,n1,n2,l1,l2,a_Array)
local G;
local count_a;
local m;
local a;
G:= Array(1..nodes,1..nodes);
for m from 1 by 1 to nodes
do
for count_a from 1 by 1 to nodes
do
a:= a_array[m,count_a];
#G[m,count_a]:= (1/omega_Two_squares(n1,a,m,l2))*((n1^2)*(int(int((Phi_left(n1,n2,l1,l2,a,m))*g,x=-l1..0,numeric=true),y=0..l2,numeric=true))) + ((n2^2)*(int(int((Phi_right(n1,n2,l1,l2,a,m))*g,x=0..l1,numeric=true),y=0..l2,numeric=true)));
G[m,count_a]:= (1/omega_Two_squares(n1,a,m,l2)) * int(int((n_of_x)^2 * PhiFunc(n1,n2,l1,l2,a,m) * g,x=-l1..l1,numeric=true),y=0..l1,numeric=true);
od;
od;
RETURN(G);
end proc:
WaveMaker_Two_squares:=proc(F,G,nodes,l1,l2,a_array,n1,n2)
local a;
local a_count;
local m;
local totalWave;
totalWave := 0;
for m from 1 by 1 to nodes
do
for a_count from 1 by 1 to nodes
do
a:= a_array[m,a_count];
totalWave := totalWave + (F[m,a_count] * wave_Two_squares(n1,n2,l1,l2,a,m)) + (G[m,a_count] * wave_Two_squares2(n1,n2,l1,l2,a,m));
end do;
end do;
RETURN(totalWave);
end proc:
AnimationMaker_Two_squares:= proc(Sigma,Vx,Vy,nodes,l1,l2,Dimensional,n1,n2,a_array,time)
local f;
local g;
local F;
local G;
local wavee;
f:= FMaker(x,y,Vx,Vy,Sigma); g:= GMaker(f,Vx,Vy,Sigma);
F:=Two_squares_CoeffFinder(f,nodes,n1,n2,l1,l2,a_array);
G:=Two_squares_CoeffFinder2(g,nodes,n1,n2,l1,l2,a_array); wavee:= WaveMaker_Two_squares(F,G,nodes,l1,l2,a_array,n1,n2);
if(Dimensional)
then
(RETURN(animate(plot3d,[wavee,x=-l1..l1,y=0..l2],t=0..time,axes=frame,frames=150, shading = zhue,scaling=constrained)));
else
(RETURN(animate(densityplot,[wavee,x=-l1..l1,y=0..l2],t=0..time,frames=150,scaling=constrained))); end if;
end proc:
2 Squares with Boundary Output
#AnimationMaker_Two_squares(10,1,0,10,Pi,Pi,false,1,2,a_array,10);
#Check initial function against f
(*
Sigma:=10;
Vx:=1;
Vy:=2;
l1:=Pi;
l2:=Pi;
n1:=1;
n2:=2;
nodes:=10;
f:= FMaker(x,y,Vx,Vy,Sigma);
g:= GMaker(f,Vx,Vy,Sigma);
F:=Two_squares_CoeffFinder(f,nodes,n1,n2,l1,l2,a_array);
G:=Two_squares_CoeffFinder2(g,nodes,n1,n2,l1,l2,a_array);
wavee:=WaveMaker_Two_squares(F,G,nodes,l1,l2,a_array,n1,n2):
plot3d(f,x=-Pi..Pi,y=0..Pi,scaling=constrained);
t:=0;
plot3d(wavee,x=-l1..l1,y=0..l1,scaling=constrained);
*)
#Check phi normalisation
(*
Vx:=1;
Vy:=2;
l1:=Pi;
l2:=Pi;
n1:=1;
n2:=2;
m:=2;
a:=1;
n_of_x:= piecewise(x<0,n1,x>0,n2);
PhiFunc:= Phi_left(n1,n2,l1,l2,a,m) *HeavyStepLeft + Phi_right(n1,n2,l1,l2,a,m) *HeavyStepRight;
int((n_of_x)^2 *(PhiFunc)^2,x=-l1..l1,y=0..l2);
*)
#Check WaveMaker
Vx:=1;
Vy:=0;
Sigma:=10;
l1:=Pi;
l2:=Pi;
n1:=1;
n2:=2;
nodes:=10;
f:= FMaker(x,y,Vx,Vy,Sigma);
g:= GMaker(f,Vx,Vy,Sigma);
F:=Two_squares_CoeffFinder(f,nodes,n1,n2,l1,l2,a_array);
G:=Two_squares_CoeffFinder2(g,nodes,n1,n2,l1,l2,a_array);
wavee:= WaveMaker_Two_squares(F,G,nodes,l1,l2,a_array,n1,n2);
animate(densityplot,[wavee,x=-l1..l1,y=0..l2],t=0..10,frames=150,scaling=constrained);
