exercise 6.12
exercise 6.12 notes:
solving PDE E''/E=-k:
1.- muPAD shows k in y''=-k*y is proportional to frequency
2.- Basic MATLAB example showing how to solve y''+abs(y)=0 y(0)=0 y(4)=-2
3.- MATLAB basic example, how to solve differential equation with parameter
4.- link to Reference [McKibben] Differential Equations with MATLAB, M McKibben M Webster, pg186.
5.- Mathworks help: Partial Differential Equations, PDE solver intro, link.
Where the solutions manual suggests to start with a function with shape E(x,y,z)=X(x)*Y(y)*Z(z) since the suggested
X Y Z are actually fields, not coordinates, let's instead use the following terms E(x,y,z)=Ex(x)*Ey(y)*Ez(z) .
SinceTM modes are assumed, both and field components of E should be null.
I tried bvp4c with syntax sol = bvp4c(odefun,bcfun,solinit) but I could not find a way to feed the function bcfun in turn used to define such input field of bvp4c with boundary conditions as parameters, and not as in-function constants defined figures only, as shown in the bvp4c intro example, that then (the boundary conditions) have to be manually modified every time a boundary condition changes. How to use bvp4c basic example shown in [2].
MuPAD shows k in y''=-k*y proportional to frequency, but MuPAD doesn't supply the general solution expected from a symbolic ODE solver, details of steps done to reach such conclusion of the limited usefulness of MuPAD in [1].
From [4]: X''(x)/X(x)=k has general solution C1*cos(n*x)+C2*sin(n*x) where k=-n^2
Actually the 3 ODEs that E ends up broken into when applying the wave equation are
Ex''/E='kx^2
Ex''/E='kx^2
Ex''/E='kx^2
where Ex is the 'phasor' of the E(x,y,z,t). The PDE that includes the time domain is already solved in [5].
To cover the 'gap' between the definition of the cavity and the assumptions
-
E(x,y,z)=X(x)*Y(y)*Z(z)
-
TM only, and
-
Ex(x)=A*cos(kx*x)+B*sin(kx*x)
-
Ex(0)=Ex(a)=0 then Ex=A*cos(kx*x), perfect conductor
-
filler dielectric ignored
all assumed in the solutions manual, and, once those 3 separate single variable PDEs solved, to get k from kx ky kz
k^2=kx^2+ky^2+kz^2 = (m*pi/a)^2+(n*pi/b)^2+(p*pi/d)^2
I recommend to use the Partial Differential Equation Toolbox. No MATLAB separate script supplied for this answer.
1.- muPAD shows k in y''=-k*y is proportional to frequency
No way to have MuPAD supplying the general solution A*cos(n*x)+B*sin(n*x), yet, although MuPAD plot does not refine the plot step, yet when frequency increases, alias corrupts the output plot, it's shown that as expected, with increasing k (wave number), decreasing wavelength, the frequency increases. MuPAD uses red font colour for MuPAD code, MuPAD's marketing manager hasn't thought it through (!?) or the logic effect after MuPAD purchased by Mathworks.
o:=ode(y''(x)=-n^2*x,y(x))
solve(o)
o2:=ode(y''(x)=-k*x,y(x))
solve(o2)
o3:=ode({y''(x)=-k*x,y(0)=0,y(a)=0},y(x))
solve(o3)
ode::solve(y''(x)=-k*x,y(x))
solve(o2,Type=ExactSecondOrder)
ode::solve(y''(x)=-k*x,y(x),Type=ExactSecondOrder)
o4:=ode({y''(x)=-k*y(x),y(0)=0,y(a)=0},y(x))
solve(o4)
solve(o4,Type=ExactSecondOrder)
f1=solve(o4)
Warning: Cannot detect the exact second-order ODE. [ode::exact_second_order]
the attempts to obtain a plot with unknown a do not render any graph:
f:=(t,Y)->[Y[2],-a*Y[1]]:Y0:=[1,0]:
G1:=(t,Y)->[t,Y[1]]:
G2:=(t,Y)->[t,Y[1]^8/8+Y[2]^2/2]:
G2:=(t,Y)->[t,Y[1]^2/2+Y[2]^2/2]:
G3:=(t,Y)->[t,Y[1]^8/8+Y[2]^2/2]:
p := plot::Ode2d(f, [i/2 $ i = 0..40], Y0,
[G1, Style = Lines, Color = RGB::Red],
[G1, Style = Points, Color = RGB::Black],
[G2, Style = Lines, Color = RGB::Blue],
[G2, Style = Points, Color = RGB::Black],
[G3, Style = Lines, Color = RGB::Green],
[G3, Style = Points, Color = RGB::Black],
PointSize = 1.5*unit::mm,
LineWidth = 0.2*unit::mm
):
now assign k=1
f := (t, Y) -> [Y[2], - Y[1]^7]:
Y0 := [1, 0]:
G1 := (t, Y) -> [t, Y[1]]:
G2 := (t, Y) -> [t, Y[1]^2/2 + Y[2]^2/2]:
G3 := (t, Y) -> [t, Y[1]^8/8 + Y[2]^2/2]:
p := plot::Ode2d(f, [i/2 $ i = 0..40], Y0,
[G1, Style = Lines, Color = RGB::Red],
[G1, Style = Points, Color = RGB::Black],
[G2, Style = Lines, Color = RGB::Blue],
[G2, Style = Points, Color = RGB::Black],
[G3, Style = Lines, Color = RGB::Green],
[G3, Style = Points, Color = RGB::Black],
PointSize = 1.5*unit::mm,
LineWidth = 0.2*unit::mm
):
plot(p):
with k=0.5
f3 := (t, Y) -> [Y[2], - 0.5*Y[1]]:
Y0 := [1, 0]:
G1 := (t, Y) -> [t, Y[1]]:
G2 := (t, Y) -> [t, Y[1]^2/2 + Y[2]^2/2]:
G3 := (t, Y) -> [t, Y[1]^8/8 + Y[2]^2/2]:
p3 := plot::Ode2d(f3, [i/2 $ i = 0..40], Y0,
[G1, Style = Lines, Color = RGB::Red],
[G1, Style = Points, Color = RGB::Black],
[G2, Style = Lines, Color = RGB::Blue],
[G2, Style = Points, Color = RGB::Black],
[G3, Style = Lines, Color = RGB::Green],
[G3, Style = Points, Color = RGB::Black],
PointSize = 1.5*unit::mm,
LineWidth = 0.2*unit::mm
):
plot(p3):
f4 := (t, Y) -> [Y[2], - 0.9*Y[1]]:
Y0 := [1, 0]:
G1 := (t, Y) -> [t, Y[1]]:
G2 := (t, Y) -> [t, Y[1]^2/2 + Y[2]^2/2]:
G3 := (t, Y) -> [t, Y[1]^8/8 + Y[2]^2/2]:
p4 := plot::Ode2d(f4, [i/2 $ i = 0..40], Y0,
[G1, Style = Lines, Color = RGB::Red],
[G1, Style = Points, Color = RGB::Black],
[G2, Style = Lines, Color = RGB::Blue],
[G2, Style = Points, Color = RGB::Black],
[G3, Style = Lines, Color = RGB::Green],
[G3, Style = Points, Color = RGB::Black],
PointSize = 1.5*unit::mm,
LineWidth = 0.2*unit::mm
):
plot(p4):
f5 := (t, Y) -> [Y[2], - 1.5*Y[1]]:
Y0 := [1, 0]:
G1 := (t, Y) -> [t, Y[1]]:
G2 := (t, Y) -> [t, Y[1]^2/2 + Y[2]^2/2]:
G3 := (t, Y) -> [t, Y[1]^8/8 + Y[2]^2/2]:
p5 := plot::Ode2d(f5, [i/2 $ i = 0..40], Y0,
[G1, Style = Lines, Color = RGB::Red],
[G1, Style = Points, Color = RGB::Black],
[G2, Style = Lines, Color = RGB::Blue],
[G2, Style = Points, Color = RGB::Black],
[G3, Style = Lines, Color = RGB::Green],
[G3, Style = Points, Color = RGB::Black],
PointSize = 1.5*unit::mm,
LineWidth = 0.2*unit::mm
):
plot(p5):
f6 := (t, Y) -> [Y[2], - 16*Y[1]]:
Y0 := [1, 0]:
G1 := (t, Y) -> [t, Y[1]]:
G2 := (t, Y) -> [t, Y[1]^2/2 + Y[2]^2/2]:
G3 := (t, Y) -> [t, Y[1]^8/8 + Y[2]^2/2]:
p6 := plot::Ode2d(f6, [i/2 $ i = 0..40], Y0,
[G1, Style = Lines, Color = RGB::Red],
[G1, Style = Points, Color = RGB::Black],
[G2, Style = Lines, Color = RGB::Blue],
[G2, Style = Points, Color = RGB::Black],
[G3, Style = Lines, Color = RGB::Green],
[G3, Style = Points, Color = RGB::Black],
PointSize = 1.5*unit::mm,
LineWidth = 0.2*unit::mm
):
plot(p6):
So, although MuPAD points on the right direction, for y''=-k*y it does not supply the expected
A*cos(sqrt(k)*x)+B*sin(sqrt(k)*x) expected solution.
However, allowing k Complex then the oscillating behaviour is hinted with
2.- Basic MATLAB example showing how to solve y''+abs(y)=0 y(0)=0 y(4)=-2
Basic example showing how to solve y''+abs(y)=0 y(0)=0 y(4)=-2
function example_04
clear all;close all;clc
x01=0 % define x range
x02=4
dx=100 % define x resolution
% define y boundary conditions
y01=0 % y(x01)=y01
y02=0 % y(x02)=y02
% build initial conditions input for bvp4c
solinit1 = bvpinit(linspace(x01,x02,dx),[1 0]);
% solve
sol1 = bvp4c(@twoode,@twobc,solinit1);
% deliver
x = linspace(x01,x02,dx)
y1 = deval(sol1,x);
figure(1);plot(x,y1(1,:));grid on
xlabel('x');ylabel('y')
% 2nd solution
solinit2 = bvpinit(linspace(x01,x02,dx),[-1 0]);
sol2 = bvp4c(@twoode,@twobc,solinit2);
x = linspace(x01,x02,dx);
y2 = deval(sol2,x);
figure(2);plot(x,y2(1,:));grid on
xlabel('x');ylabel('y')
% support functions
function res = twobc(ya,yb)
res = [ ya(1); yb(1)+2 ];
function dydx = twoode(x,y)
dydx = [ y(2); -abs(y(1)) ];
3.- MATLAB basic example, how to solve differential equation with parameter
equation: y''+(lambda-2*q*cos(2*x))*y=0
boundary: y(0)=1 y'(0)=0 y'(pi)=0
either place the following along with the 3 support functions in a single file with function header containing all support functions, or define the support functions in same folder as the following lines.
lambda = 15;
solinit = bvpinit(linspace(0,pi,10),@mat4init,lambda);
sol = bvp4c(@mat4ode,@mat4bc,solinit);
fprintf('The fourth eigenvalue is approximately %7.3f.\n', sol.parameters)
xint = linspace(0,pi);
Sxint = deval(sol,xint);
plot(xint,Sxint(1,:))
axis([0 pi -1 1.1])
title('Eigenfunction of Mathieu''s equation.')
xlabel('x')
ylabel('solution y')
function dydx = mat4ode(x,y,lambda)
q = 5;
dydx = [ y(2)
-(lambda - 2*q*cos(2*x))*y(1) ]
function yinit = mat4init(x)
yinit = [ cos(4*x)
-4*sin(4*x) ];
function res = mat4bc(ya,yb,lambda)
res = [ ya(2)
yb(2)
ya(1)-1 ];
4.- link to Reference [McKibben] Differential Equations with MATLAB,
Mark McKibben Micah Webster, pg186.
5.- Link to Mathworks help: Partial Differential Equations, PDE solver
https://uk.mathworks.com/help/matlab/math/partial-differential-equations.html?s_tid=srchtitle