top of page

exercise 6.12

 exercise 6.12 notes:

Small tag OK.jpg
pozar_06_exercise_12_question.jpg

solving PDE E''/E=-k:

Small tag OK.jpg

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.

001_1.jpg
001_2.jpg
001_3.jpg

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)

002_1.jpg
002_2.jpg
002_3.jpg
002_4.jpg
002_5.jpg
002_6.jpg

Warning: Cannot detect the exact second-order ODE. [ode::exact_second_order]

002_7.jpg

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

003_1.jpg

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

​

​

003_2.jpg
003_3.jpg
003_4.jpg
003_5.jpg
003_6.jpg

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

002_7.jpg

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

001.jpg
002.jpg

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.

​

https://www.amazon.co.uk/Differential-Equations-MATLAB-Exploration-Applications-ebook/dp/B00NFODJAQ/ref=sr_1_fkmrnull_1?keywords=mckibben+matlab+differential+equations&qid=1558098889&s=gateway&sr=8-1-fkmrnull

​

​

5.- Link to Mathworks help: Partial Differential Equations, PDE solver

​

https://uk.mathworks.com/help/matlab/math/partial-differential-equations.html?s_tid=srchtitle

bottom of page