Review: Numerical Methods for Differential Equations

ECE Generic and 2331

Han Q. Le(c)

Segment 3 - Practice and examples - ODE

1. Basic use of NDSolve

Example how to use NDSolve

s=NDSolve[{y'[x]==y[x]Cos[x+y[x]],y[0]==1},y,{x,0,30}]

Plot[Evaluate[y[x]/.s],{x,0,30},PlotRange->All]

2. Stable pendulum

2.1 Single pendulum

Go to webpage to download the full app. Here shows only partial results to make the code simple to understand.

code

Block[{T,sol,xarr},
Manipulate[
If[TrueQ[RUNpendum1== ],
Manipulate[
T=(4/(2 π (Sqrt[2]))) NIntegrate[(1./(Sqrt[Cos[u]-Cos[θinit]])),{u,0.,θinit}] ;
(* T=Chop[(4/(2 π (Sqrt[2])))((2.` EllipticF[((θinit)/2),((Csc[((θinit)/2)])^2)])/(Sqrt[1-Cos[θinit]])) ]; *)

sol=NDSolve[{x''[t]+(2 π)^2 Sin[x[t]]==0,x[0]==θinit,x'[0]==vinit},x,{t,T}
,Method-> {"ExplicitRungeKutta" }];
xarr=Table[Evaluate[x[t]/.First[sol]],{t,0,199T/200.,T/200.}];
Manipulate[
Graphics[{Hue[0.5,1.,1.],Thickness[0.01],Line[{{0.,0.},{Sin[xarr[[i]] ],-Cos[xarr[[i]]  ]}}]
,White,PointSize[0.02],Point[{0.,0.}],
Red,Disk[{Sin[xarr[[i]] ],-Cos[xarr[[i]]  ]},0.1]
,Inset[NumberForm[(i-1)*T/100.,2],ImageScaled[{0.9,0.9}]
,BaseStyle-> {24, Yellow}]}
,PlotRange-> {{-1.25,1.25},{-1.25,1.25}}
,Background-> Black]
,{i,1,Length[xarr],1,ControlType-> Animator, DefaultDuration-> 2T}]
,{{θinit,Pi/4.},0.05,0.95Pi},{{vinit,0.},-10,10.}]
,Button[,RUNpendum1= ]]
,
Control[{{RUNpendum1,,},{,},ControlType->Toggler}] ]
]

code

Block[{deqns,aeqns,ics,sol1},
Manipulate[
If[TrueQ[RUNpendum2== ],
Manipulate[
sol1=NDSolve[{x''[t]+(((2 π))^2) Sin[x[t]]==0,x'[0]==vinit,x[0]==θinit},x,{t,0,5}];
Plot[Evaluate[{x[t],x'[t]}/.sol1],{t,0,5}]
,{{θinit,Pi/2.},0.05,0.95Pi},{{vinit,0.},-10,10.}
,Initialization:>{
deqns={x''[t]==(((2 π))^2) Sin[x[t]]};
aeqns={((x[t])^2)+((y[t])^2)==(1^2)};
ics={x[0]==1,y'[0]==1};}]
,Button[,RUNpendum2= ]]
,
Control[{{RUNpendum2,,},{,},ControlType->Toggler}] ]
]

2.2 Double pendulum

code

Block[{deqnsDP,aeqnsDP,icsDP,solDP},
Manipulate[
If[TrueQ[RUNdoublependul== ],
Manipulate[
icsDP={x1[0]==x1init,x2[0]==x2init,y1'[0]==y1pinit,y2'[0]==y2pinit};
solDP=Quiet@NDSolve[{deqnsDP,aeqnsDP,icsDP},{x1,y1,x2,y2,λ1,λ2},{t,0,15},Method->{"IndexReduction"->{True,"ConstraintMethod"->"Projection"}}] ;
Animate[Graphics[{{Thickness[0.02],Hue[0.15,1,1],Line[{{0,0},{x1[t],y1[t]}}]},{Thickness[0.02],Hue[0.5,1,1],Line[{{x1[t],y1[t]},{x2[t],y2[t]}}]}}/.solDP,PlotRange->{{-2.1,2.1},{-2.1,2.1}},Frame -> True
,Background->Black],{t,0,15},SaveDefinitions->True,AnimationRunning->False]
,{{x1init,1},0.,1},{{x2init,1},0.,2}
,{{y1pinit,0},0.,0.5},{{y2pinit,0.},0.,0.5}]

,Button[,RUNdoublependul= ]]
,
Control[{{RUNdoublependul,,},{,},ControlType->Toggler}] 
,Initialization:>{
deqnsDP={x1''[t]==λ1[t] x1[t]-λ2[t](x2[t]-x1[t]), y1''[t]==λ1[t] y1[t]-λ2[t](y2[t]-y1[t])-9.81,
x2''[t]==λ2[t](x2[t]-x1[t]),y2''[t]==λ2[t] (y2[t]-y1[t])- 9.81};
aeqnsDP={((x1[t])^2)+((y1[t])^2)==1,(((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2) ==1};

}]
]

2.3 Chaos concept

Class discussion

Additional demo

https://www.youtube.com/watch?v=QXf95_EKS6E

3. Inverted pendulum

3.1 Introduction

Concept of stability. See below.

We now consider a well-studied nonlinear problem of stabilization. But here, we will only investigate the linearized case. Consider an inverted pendulum problem below.

There are only 2 degrees of freedom: x and θ as shown. But the pendulum moves relative to the cart. Newton 3rd law stipulates that the force (Subscript[F, x]), (Subscript[F, y]) the cart exerts on the pendulum has to be equal and opposite to that of the pendulum on the cart. Hence, the equation of motion is:
For the cart center of mass:       M (Overscript[x, ..])= u -(Subscript[F, x])                       (4.1.12)
For the pendulum center of mass:         m ((d^2)/(d(t^2))) (x -l Sin[θ])= (Subscript[F, x])      (4.1.13)
                                                              m ((d^2)/(d(t^2))) (l Cos[θ])= -m g +(Subscript[F, y])      (4.1.14)
For the pendulum rotation motion:        J (Overscript[θ, ..]) =(Subscript[F, x]) l Cos[θ]+(Subscript[F, y]) l Sin[θ]      (4.1.15)
Note: the polarity of θ is counter-clockwise (right hand rule).

We have 4 Eqs. for 4 unknowns. But we are not interested in (Subscript[F, x]) and (Subscript[F, y]), hence we eliminate them:
(Combine 4.1.13 into 4.1.12):
                        M (Overscript[x, ..])= u -(Subscript[F, x])=u-m ((d^2)/(d(t^2))) (x -l Sin[θ])      (4.1.16)
 and combine both (4.1.13) and (4.1.15) in (4.1.16):
                                      J (Overscript[θ, ..]) =(Subscript[F, x]) l Cos[θ]+(Subscript[F, y]) l Sin[θ]                              (4.1.15) 
                                      J (Overscript[θ, ..]) =m ((d^2)/(d(t^2))) (x-l Sin[θ]) l Cos[θ]
              +(m g+m ((d^2)/(d(t^2))) (l Cos[θ]) ) l Sin[θ]     (4.1.17)

The final eqs are:              M (Overscript[x, ..])= u-m ((d^2)/(d(t^2))) (x -l Sin[θ])      (4.1.16)
                           J (Overscript[θ, ..]) =m ((d^2)/(d(t^2))) (x-l Sin[θ]) l Cos[θ] +(m g+m ((d^2)/(d(t^2))) (l Cos[θ]) ) l Sin[θ]     (4.1.17)

Or:                                           (Overscript[x, ..])= (u/M)-(m/M) ((d^2)/(d(t^2))) (x -l Sin[θ])      (4.1.16)
                                                (Overscript[θ, ..]) =((m (l^2))/J) ((d^2)/(d(t^2))) ((x/l)-Sin[θ]) Cos[θ] +((m (l^2))/J)((g/l)+ ((d^2)/(d(t^2))) Cos[θ] ) Sin[θ]

Or:                                           ((Overscript[x, ..])/l)= ((s[t])/(τ^2))-μ ((d^2)/(d(t^2))) ((x/l) - Sin[θ])      (4.1.16)
                                                (Overscript[θ, ..]) =ν ( ((d^2)/(d(t^2))) ((x/l)-Sin[θ]) Cos[θ] +(((Subscript[ω, 0])^2)+ ((d^2)/(d(t^2))) Cos[θ] ) Sin[θ])    
                                               let J=((m (l^2))/12) then  ν=3.

3.2 Linearization

Now, we have 2 Eqs (4.1.16), (4.1.17) with 2 unknowns x and θ, but they are transcendental and very difficult to solve. We linearize the equation by assuming that θ is a very small angle, barely deviated from the balance vertical position. Then, Cos[θ] ~ 1 and Sin[θ] ~ θ. The Eqs become:
                              M (Overscript[x, ..])=u -m ((d^2)/(d(t^2))) (x -l θ)
        =u-m (Overscript[x, ..])+m l(Overscript[θ, ..])     
Or:                          (M+m) (Overscript[x, ..]) -m l(Overscript[θ, ..])=u                               (4.1.18) 
and:                     J (Overscript[θ, ..]) =m ((d^2)/(d(t^2))) (x -l θ ) l+m g l θ
         = m l (Overscript[x, ..])-m (l^2) (Overscript[θ, ..])+m g l θ      
Or:                   m l (Overscript[x, ..]) - (J+m (l^2))(Overscript[θ, ..]) +m g l θ =0            (4.1.19)

We apply LT and rewrite (4.1.18) and (4.1.19) in matrix form:
                     (({{(M+m) (s^2), -m l (s^2)}, {m l (s^2), -(J+m (l^2)) (s^2)+m g l}})) . (({{X}, {Θ}}))==(({{U}, {0}}))      (4.1.20)

Solve[ (({{(M+m) (s^2), -m l (s^2)}, {m l (s^2), -(J+m (l^2)) (s^2)+m g l}})) . (({{X}, {Θ}}))==(({{U}, {0}})) ,{X,Θ}]

FullSimplify[Inverse[(({{(M+m) (s^2), -m l (s^2)}, {m l (s^2), -(J+m (l^2)) (s^2)+m g l}}))] ]

The solution is:
(({{X}, {Θ}}))=(({{(((m (l^2)+J) (s^2)-g l m)/((m M (l^2)+J (m+M)) (s^4)-g l m (m+M) (s^2))), ((l m)/(g l m (m+M)-(m M (l^2)+J (m+M)) (s^2)))}, {((l m)/((m M (l^2)+J (m+M)) (s^2)-g l m (m+M))), ((m+M)/(g l m (m+M)-(m M (l^2)+J (m+M)) (s^2)))}})) . (({{1}, {0}})) U    (4.1.21)

Simplify[Det[(({{(M+m) (s^2), -m l (s^2)}, {m l (s^2), -(J+m (l^2)) (s^2)+m g l}}))]]

Define:           Q[s]=-(s^2) ((s^2) (J (m+M)+(l^2) m M)-g l m (m+M))       (4.1.22)
We see that Q[s] has just a double root: s=0 and
                                   (s^2)=((g l m (m+M))/(J (m+M)+(l^2) m M)) ;                       (4.1.23)
We define:   (Sqrt[((g l m (m+M))/(J (m+M)+(l^2) m M))])≡α                                     (4.1.24)

Then:
(({{X}, {Θ}}))=(({{(((m (l^2)+J) (s^2)-g l m)/((m M (l^2)+J (m+M)) (s^4)-g l m (m+M) (s^2)))}, {((l m)/((m M (l^2)+J (m+M)) (s^2)-g l m (m+M)))}}))  U  
         =(({{(((m (l^2)+J) )/((m M (l^2)+J (m+M))))  (((s^2)-(β^2))/((s^2) ( (s^2)-(α^2))))}, {((l m)/((m M (l^2)+J (m+M))))(1/(( (s^2)-(α^2)))) }})) U                                (4.1.25a)                
 where we define: (( g l m)/((m (l^2)+J) ))≡ (β^2)                                          (4.1.25b)
      =(({{((g l  m (m+M) )/((m M (l^2)+J (m+M))))(((m (l^2)+J))/(g l m (m+M) )) (((s^2)-(β^2))/((s^2) ( (s^2)-(α^2))))}, {((g l m (m+M))/((m M (l^2)+J (m+M)))) ((l m)/(g l m (m+M))) (1/(( (s^2)-(α^2)))) }})) U  
(({{X}, {Θ}}))=(({{(1/((m+M)))((α^2)/( ((β )^2))) (((s^2)-(β^2))/((s^2) ( (s^2)-(α^2))))}, {(1/( g (m+M))) ((α^2)/(( (s^2)-(α^2)))) }})) U                                          (4.1.26)

Thus, the TF are:
                      (({{(X/U)}, {(Θ/U)}}))=(({{(1/((m+M)))((α^2)/( ((β )^2))) (((s^2)-(β^2))/((s^2) ( (s^2)-(α^2))))}, {(1/( g (m+M))) ((α^2)/(( (s^2)-(α^2)))) }}))                     (4.1.27)

3.2.1 Step response

If the pendulum is perfectly balance at vertical position, and we knock the cart with a force. What happens?

Let’s find out:

FullSimplify[InverseLaplaceTransform[ (1/s) ((α^2)/(( (s^2)-(α^2)))),s,t]]

Hence, the solution for θ is:
θ[t]=(u/( g (m+M))) (cosh(α t)-1) 
This is valid for only small θ, hence, we consider an expansion for small α t:
θ[t]~(u/( g (m+M))) (1+((((α t))^2)/2)+((((α t))^4)/(4!))+...-1) =((u (α^2) )/( 2 g (m+M)))((t^2) +((α^2)/12)(t^4)+...) 
What this means is that θ rapidly grows and it doesn’t return to the original value of 0. This is an unstable case.

To find the solution for x.

FullSimplify[InverseLaplaceTransform[(((α^2)/( ((β )^2))))  (1/s)(((s^2)-(β^2))/((s^2) ( (s^2)-(α^2)))),s,t]]

((t^2)/2)+FullSimplify[(1/(α^2))-(1/(β^2))+((1/(β^2))-(1/(α^2))) cosh(α t)]

FullSimplify[(((α^2)-(β^2))/((α^2) (β^2)))  /. α->(Sqrt[ ((g l m (m+M))/(J (m+M)+(l^2) m M))]) /. β-> (Sqrt[((g l m)/((m (l^2)+J) ))])]

The solution for x is:
x[t]=(1/2)(u/(M+m))(t^2) +((l m u)/(g (((m+M))^2))) (cosh(α t)-1)
Again, we see that the term (1/2)(u/(M+m))(t^2) represent the effect of force u with acceleration (u/(M+m)). But the term ((l m u)/(g (((m+M))^2))) (cosh(α t)-1) represents the effect of another force on it: this is the kick by the pendulum as it swings downward.

Manipulate[
Show[Graphics[{Red,Rotate[Rectangle[{xIP[t,{α1,β1}]-0.05,0.67}
,{xIP[t,{α1,β1}]+0.05,2.67}],θIP[t,α1],{xIP[t,{α1,β1}],0.67} ]
}] 
,Graphics[{Blue,Rectangle[{xIP[t,{α1,β1}]-0.5,0},{xIP[t,{α1,β1}]+0.5,0.67}]}] 
, PlotRange->{ {-1.25,1.25},{0,3}}]
,{t,0,1.5}
,Initialization:>{
 α1= 0.5 ; β1= 0.4 ; 
xIP[t_,{α_,β_}]:=((t^2)/2)+(((α-β) (α+β) )/((α^2) (β^2)))(Cosh[α t]-1) ;
θIP[t_,α_]:=Cosh[α t]-1. ;l=1. ;

}]

Let’s take a look at the roots of the TFs:
                      (({{(X/U)}, {(Θ/U)}}))=(({{(1/((m+M)))((α^2)/( ((β )^2))) (((s^2)-(β^2))/((s^2) ( (s^2)-(α^2))))}, {(1/( g (m+M))) ((α^2)/(( (s^2)-(α^2)))) }}))                     (4.1.27)

For simplicity, let’s consider the case of uniform diameter and uniform density rod. Moment of inertia is:

(m/((π (ρ^2)2 l)))(∫_(-l)^l) π (ρ^2)(x^2) \[DifferentialD]x

α=(Sqrt[((g l m (m+M))/(J (m+M)+(l^2) m M))]) /. J->(1/3) m (l^2)

FullSimplify[(Sqrt[((g l m (m+M))/((l^2) m M+(1/3) (l^2) m (m+M)))])]

Manipulate[ g=1.;
α=(Sqrt[((3 g (m+M))/(l (m+13 M)))]); 
CPPlot[{-α,0,α},2,{0.01,0.05,0.05},,]
,{{m,0.5,"pendulum mass"},0.1,5},{{M,5,"carrier mass"},1,20},{{l,1,"pendulum length"},0.2,5}]

Notice the root in the real positive half plane

3.2.2 Step response

What if the force is harmonic? From (4.1.26)

(({{X}, {Θ}}))=(({{(1/((m+M)))((α^2)/( ((β )^2))) (((s^2)-(β^2))/((s^2) ( (s^2)-(α^2))))}, {(1/( g (m+M))) ((α^2)/(( (s^2)-(α^2)))) }})) U              (4.1.26)

We can solve for the harmonic case by simply substitute i ω for s

(({{(1/((m+M)))((α^2)/( ((β )^2))) (((s^2)-(β^2))/((s^2) ( (s^2)-(α^2))))}, {(1/( g (m+M))) ((α^2)/(( (s^2)-(α^2)))) }}))  /. s-> I ω

Hence, the solutions are:
x[t]=-(((α^2) ((β^2)+(ω^2)))/((β^2) (ω^2) ((α^2)+(ω^2))))((Subscript[u, 0])/(m+M))Sin[ω t]            (4.1.27a)
θ[t]=-((α^2)/(((α^2)+(ω^2))))((Subscript[u, 0])/(g (m+M))) Sin[ω t]                (4.1.27b)
As long as   ((α^2)/(((α^2)+(ω^2))))((Subscript[u, 0])/(g (m+M))) is small, the solution is approximately valid (although true stability proof requires the full consideration of the transcendental equations). In order words, active force can be applied such that θ is within bound (stable), something that we see people and animals can do.

3.3 Numerical calculation

See HW.

Additional info and demo

http://www.youtube.com/watch?v=MWJHcI7UcuE
http://www.youtube.com/watch?v=AZhQt7HOSWo
http://www.youtube.com/watch?v=a4c7AwHFkT8
http://www.youtube.com/watch?v=cyN-CRNrb3E
https://youtu.be/15DIidigArA
http://www.youtube.com/watch?v=bcr3cG2AKBs
http://www.youtube.com/watch?v=h-oZoNC_SIo
http://www.youtube.com/watch?v=LTyGdRG8y9Q

4. Demo from Mathematica: bouncing ball (energy attenuation)

4.1 Bouncing ball

Hybrid Differential Equation.
Simulate a bouncing ball that retains 95% of its velocity in each bounce:

NDSolve[{y''[t]==-9.81,y[0]==5,y'[0]==0,WhenEvent[y[t]==0,y'[t]->-0.95y'[t]]},y,{t,0,10}];

Plot[y[t]/.%,{t,0,10}]

Model a ball bouncing down steps:

ballsteps=NDSolve[{x''[t]==0,y''[t]==-9.8,y[0]==6,y'[0]==0,x[0]==0,x'[0]==1,a[0]==5,WhenEvent[Mod[x[t],1]==0,If[a[t]>0,a[t]->a[t]-1,"RemoveEvent"]],WhenEvent[y[t]==a[t],{x'[t],y'[t]}->.9{x'[t], -y'[t]}]},{x,y,a},{t,0,15},DiscreteVariables->{a}];

Show[ParametricPlot[Evaluate[{{x[t],y[t]},{x[t],a[t]}} /. ballsteps],{t,0,15}], Plot[{0,Floor[6-x]},{x,-1,15}, Filling->{2->0}, Exclusions->None],Frame->{{True,False},{True,False}},FrameLabel->{"x","y"}]

Each time a linear oscillator solution crosses the negative x axis, reflect it across the y axis:

de[δ_]:={(Subscript[x, 1])'[t]==(Subscript[x, 2])[t],(Subscript[x, 2])'[t]==-(Subscript[x, 1])[t]+2δ (Subscript[x, 2])[t]+1};
ic={(Subscript[x, 1])[0]==.5,(Subscript[x, 2])[0]==.5};

{sol,{points}}=NDSolve[{de[0.1],ic,WhenEvent[And[(Subscript[x, 2])[t]==0,(Subscript[x, 1])[t]<0],{(Subscript[x, 1])[t]->-(Subscript[x, 1])[t],Sow[-(Subscript[x, 1])[t]]}]},{(Subscript[x, 1]),(Subscript[x, 2])},{t,0,1000}]//Reap;

The solution of this reset oscillator exhibits chaotic behavior:

Map[ParametricPlot[{(Subscript[x, 1])[t], (Subscript[x, 2])[t]}/.sol, {t, 0, #}, PlotPoints->500,PlotRange->{{-1,3},{-2,2}}]&,{10,100,400}]

Plot the solution on the negative x axis with a histogram of the reflection points:

Grid[{{Histogram[points,{-1,0, .05}]},{ParametricPlot[{(Subscript[x, 1])[t], (Subscript[x, 2])[t]}/.sol, {t, 0, 1000}, PlotRange->{{-1,0},{-0.3,0}},PlotPoints->2000,Axes->None]}}]

5. Non-linear response - Op amp oscillator

This problem is in your HW. Please go to this web site for details:
http://courses.egr.uh.edu/ECE/ECE3340/Class%20 Notes2100/Lab_5/ECE2100_lab%20 V_PID.htm
(scroll down to see the above circuit)

You can download the app to check for your answer. You are NOT required to do animation, only to solve the circuit and plot the time response for a few periods.

6. Non-linear response - amplifier

The input voltage ve varies sinusoidally:

ve[t_]:=(4/10) Sin[200 Pi t];
params={vb->6,r0->1000,r1->9000,r2->9000,r3->9000,r4->9000,r5->9000,α->99/100,β->10^-6,c1->10^-6,c2->2 10^-6,c3->3 10^-6};

The transistor dispatches the voltage in a nonlinear way, depending on (Subscript[v, 2])-(Subscript[v, 3]):

v23[t_]=β (Exp[(v2[t]-v3[t])/.026]-1);

Use Ohm's law and Kirchoff's law to determine the governing equations for each node:

node1=c1 (v2'[t]-v1'[t])==v1[t]/r0-ve[t]/r0;
node2=c1(v1'[t]-v2'[t])==(1-α)v23[t]+v2 [t]/r1+v2[t]/r2-vb/r2;
node3=c2 v3'[t]==v23[t]-v3[t]/r3;
node4=c3(v4'[t]-v5'[t])==vb/r4-v4[t]/r4-α v23[t];
node5=c3(v4'[t]-v5'[t])==v5[t]/r5;
ics={v1[0]==0,v2[0]==vb/2,v3[0]==vb/2,v4[0]==vb,v5[0]==0};

sol=.;

sol=NDSolve[{node1,node2,node3,node4,node5,ics}/.params,{v1,v2,v3,v4,v5},{t,0,0.5},Method->{"EquationSimplification"->"Residual"}, MaxSteps->200000] ;

The transistor amplifies voltage (Subscript[v, 4]) relative to (Subscript[v, 3]):

Plot[Evaluate[{v3[t],v4[t]}/.sol],{t,0,0.5},PlotRange->All]

Plot[Evaluate[{v4[t]}/.sol],{t,0,0.2},PlotRange->All]

v4arr=Flatten[Table[v4[t]/.sol,{t,0.,0.5,0.001}]] ;
{v4PSD,W,{fr,v4spec}}=PSDFunc[v4arr,0.001] ;

ListLogPlot[v4PSD, PlotRange-> All, Joined-> True
,PlotStyle-> {Thick,Red}
,Filling-> {1-> {Bottom,Hue[0.05,1.,1.,0.3]}}
,Frame-> True, GridLines-> Automatic
,FrameLabel-> (Style[#,16,Black,Bold, Opacity[1.],FontFamily-> "Arial"] & /@ {"Frequency","Spectral Power Density of NRFS"})
,LabelStyle-> {16,Black,Opacity[1.],FontFamily-> "Arial"}

,ImageSize-> 600]

Segment 4 - PDE

6. FDTD of EM wave in 1-dimensional air-dielectric

FDTD

Source code for FDTD demo: pulse in air and dielectric media

Manipulate[
If[TrueQ[RUNFDTD],

Manipulate[
(* If[!TrueQ[{nold,κold}=={n1, κ}],newparm=True;{nold,κold}={n1, κ}
,]; *)

Column[{
If[!TrueQ[iready==1],
medium=Graphics[{Black,Line[{{0,0},{1600,0}}]
,Hue[0.32,((n1-1)/2.8),1.], Opacity[0.3+0.2*((n1-1)/2.8)],Rectangle[{800,-2},{1600,2}]}
,AspectRatio-> 0.2, ImageSize-> 600
,PlotRange-> {{0,1600},{-2,2}}]
,]
,If[ TrueQ[CALC] ,
n1=Max[n1,1.];
ϵ0= (n1+I κ)^2  ;
pulseT=Max[Min[τ,400],150] ;
S=Table[Quiet[(E^(-I 2. π (i-300.)/50.)) ]*If[(((i-800.))^2)/pulseT^2>8.,0.,Exp[-(((i-800.))^2)/pulseT^2]],{i,1,3000}] ;
F=ConstantArray[0.,1600] ;
H=ConstantArray[0.,1600] ;
ϵArr=Table[If[i<800,1.,ϵ0],{i,2,Length[H]}];  
 {AH,AF}=Reap[For[nt=1,nt≤5000,nt++,
F1=F ;
PrependTo[F,If[nt>1000,0.,S[[nt]] ] ] ;
F0=Drop[F,-1] ; H0=H ; 
H=H0-0.5*(F1-F0) ;
Sow[H,ah] ;
H1=Drop[H,1];H0=Drop[H,-1] ;F2=Drop[F1,-1] ;
F=F2-0.5*(H1-H0)/ϵArr ;
AppendTo[F,0.] ;
Sow[F,af]
]][[2]] ;
iready=1; CALC=False
,]

,If[TrueQ[iready==1],

Manipulate[
Show[ListPlot[{Re[AF[[i]]]+1,Re[AH[[i]]]-1}, Joined-> True, PlotRange-> {All,{-2,2}}, PlotStyle-> {{Thick,Red},{Thick,Blue},{Thick,Green},{Thick,Yellow}}
, AspectRatio-> 0.2],medium
,ImageSize-> 600],{i,1,Length[AF],1}]
,]
}]
,Row[{
Column[{Button["Medium",BaseStyle-> {14, Bold,Black,FontFamily-> "Arial"}
,Background-> Hue[0.16,0.7,1]]
,
Control[{{n1,2,"index of refraction"},1.3,3.5,ImageSize-> Medium}]
,Control[{{κ,0.,"imaginary index"},0.,0.01,ImageSize-> Medium}]
,Control[{{τ,150,"pulse duration"},150,400,ImageSize-> Medium}]
},Frame-> True
,FrameStyle->{Thick,Hue[0.,1.,1]}
,Background-> Hue[0.1,0.1,1]]
,Spacer[25]
,Grid[{{
Button["Click to run calculation",{CALC=True},BaseStyle-> {14, Bold,White,FontFamily-> "Arial"}
,Background-> Hue[0.65,1,1]]}
,{Button[Dynamic[If[nt<5001,nt,"DONE"]],BaseStyle-> {14, Bold,Black,FontFamily-> "Arial"}
,Background-> Hue[0.9,0.3,1]]}
}]
(* ,Column[{Button["generic control 2",BaseStyle-> {14, Bold,Black,FontFamily-> "Arial"}
,Background-> Hue[0.6,0.5,1]]
,
Control[{{var3,var3init,"var 3 name"},0.,1,ImageSize-> Medium}]
,Control[{{var4,var4init,"var 4 name"},0.,1.,ImageSize-> Medium}]

},Frame-> True
,FrameStyle->{Thick,Hue[0.65,1.,1]}
,Background-> Hue[0.55,0.1,1]] *)
}]



]
(* ,Manipulate[If[TrueQ[TEXTINST==on],TextCell["this is text area for notes and explanations" 
,Black,If[NumberQ[ftsz],ftsz,1]*4+14,Italic,FontFamily->"Calibri",Background->Hue[0.42,0.3,1]],""]
,Grid[{{Button["Text ->",
{If[!TrueQ[TEXTINST==on],TEXTINST=on,TEXTINST=off]}
,BaseStyle->{20,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.25,1,1]]
,Control[{{TEXTINST,on,""},{off,on},ControlType-> Toggler}]
,Button[,{If[NumberQ[ftsz], ftsz=Mod[ftsz+1,3],ftsz=1];TEXTINST=on}
,Background->  Hue[0.1,0.7,1.]]}}]
] *)
,Button[Show[,ImageSize-> 500]
,{RUNFDTD=True ;subflag=1}]]

,Column[{(* Show[,ImageSize-> {400,35}]  *)
Show[,ImageSize-> 500]
,Button["ECE generic - Demo FDTD 1-dimensional",
BaseStyle->{18,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.4,1,1]]}]
,Row[{Button["RUN STATUS/CONTROL",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUNFDTD,False,"->"},{False,True},ControlType-> Checkbox}] 
,Spacer[35]
,Button[Dynamic[If[TrueQ[RUNFDTD],,]]
,{ If[TrueQ[RUNFDTD],RUNFDTD=False,RUNFDTD=True]}]}]
,Initialization:>{
Off[General::munfl];
iready=0;nt=0;

F=ConstantArray[0.,1600] ;
H=ConstantArray[0.,1600] ;
{off,on}={Button["OFF",BaseStyle->{16, Bold,White,FontFamily-> "Arial"}
,Background->Red]
,Button["ON",BaseStyle->{16, Bold,White,FontFamily-> "Arial"}
,Background->Hue[0.36,1,0.4]]}
}
]

Demo

7. Nonlinear wave: soliton

7.1 soliton

code for soliton with analytic solution

Manipulate[Plot[Re[ESol[1.55, SolParm, z,t]],{z,-10,100}, PlotRange->{-1,1},
PlotStyle-> {Hue[0.75,0.7,1]}
, PlotPoints-> 200, ImageSize-> 700, AspectRatio->0.2],{t,0,0.5,0.0025}
,Initialization:>{
SolParm={1.45, 1.8,20. ,0.02}  ;
ESol[λ_, SolParm_List,z_, t_]:=Module[{neff,ng, zc, tc, c=299.792458} ,
{neff,ng, zc, tc}=SolParm  ;
(E^(I 2 π  (neff  z /λ - t  c/ λ))) (E^((I z)/(2 zc))) Sech[((t-((z  ng)/( c)))/tc)]   ] ;
} ]

Demo

more about soliton

http://www.youtube.com/watch?v=Ud7STKWNmQw 
http://www.youtube.com/watch?v=MADng1fqECY

http://www.youtube.com/watch?v=SAbQ4MvDqEE

http://www.youtube.com/watch?v=wEbYELtGZwI

7.2 Numerical calculation for soliton

See HW.

8. Heat diffusion (not for ECE3340)