(*   Solution of the Wave Equation using the method of separation of variables  *)

(*   Section 1 from the notes. Dirichlet Boundary Conditions   *)

(*   Consider the wave equation (Subscript[u, tt]) = (c^2) (Subscript[u, xx]), where 0<x<l and 0<t<t_0   *)
(*   with conditions u(0,t)=0   u(l,t)=0   u(x,0)=ϕ(x)  (Subscript[u, t])(x,0)=ψ(x)   *)
(*   Solve the above problem using the method of separation of variables u(x,t)=X(x)T(t) and plot the solution. *)

c:=1; l:=1;t0:=3;

(*   Excercise 4 in the notes   *)

(Subscript[f, 1])[n_,t_]:=Cos[((n*Pi*c*t)/l)]; (Subscript[f, 2])[n_,t_]:=Sin[((n*Pi*c*t)/l)]     (*  functions cos[((nπct)/l)] and sin[((nπct)/l)] solving the problem for T(t)   *)

Graph1[n_]:=Plot[{(Subscript[f, 1])[n,t],(Subscript[f, 2])[n,t]}, {t,0,t0}, PlotLabel-> "(Subscript[f, 1])(t)=cos[nπct/l]    and    (Subscript[f, 2])(t)=sin[nπct/l]"]

Graph1[1]

Graph1[2]

Graph1[3]

X[n_,x_]:=Sin[n*Pi*x/l]       (*   solution of the eigenvalue problem for X(x)   *)

Graph2[n_]:=Plot[X[n,x],{x,0,l}, PlotLabel->"(Subscript[X, n])(x)=sin[nπx/l]"]

Graph2[1]

Graph2[2]

Graph2[3]

(* Example 1. *)

ϕ[x_]:= x(1-x); ψ[x_]:=0; (*   initial conditions   *)

Plot[{ϕ[x]},{x,0,l}, PlotLabel -> "ϕ(x)" ]

A[n_]:=2/l * Integrate[ϕ[x]*Sin[n*Pi*x], {x,0,l}]; B[n_]:=2/(n*Pi*c) * Integrate[ψ[x]*Sin[n*Pi*x], {x,0,l}];  (*   coefficients (Subscript[A, n]) and (Subscript[B, n]) in the solution u(x,t)   *)

A[1]

A[2]

A[3]

A[4]

A[5]

A[6]

A[7]

(*   note that all B_n's are zero   *)

u[n_, x_, t_]:=(A[n]*(Subscript[f, 1])[n,t]+B[n]*(Subscript[f, 2])[n,t])*X[n,x]

u[1,x,t]

Plot3D[((8 Cos[π t] Sin[π x])/(π^3)), {t,0,t0}, {x,0,l}, PlotLabel -> "(Subscript[u, 1])(x,t)"]         (* (Subscript[u, 1])(x,t) *)

u[2,x,t]

u[3,x,t]

Plot3D[((8 Cos[3 π t] Sin[3 π x])/(27 (π^3))), {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 3])(x,t)"]   (* (Subscript[u, 3])(x,t) *)

u[5,x,t]

Plot3D[((8 Cos[5 π t] Sin[5 π x])/(125 (π^3))), {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 5])(x,t)"]   (* (Subscript[u, 5])(x,t) *)

u[7,x,t]

Plot3D[((8 Cos[7 π t] Sin[7 π x])/(343 (π^3))), {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 7])(x,t)"]   (* (Subscript[u, 7])(x,t) *)

Plot3D[((8 Cos[π t] Sin[π x])/(π^3))+((8 Cos[3 π t] Sin[3 π x])/(27 (π^3)))+((8 Cos[5 π t] Sin[5 π x])/(125 (π^3)))+((8 Cos[7 π t] Sin[7 π x])/(343 (π^3))), {t,0,t0}, {x,0,l}, PlotLabel->"u(x,t)=(Subscript[(∑^8), n = 1])(Subscript[u, n])(x,t) "]     (*   solution u(x,t)   *)

(*  Example 2.  *)

ϕ[x_]:= Sin[5Pi*x] + 2 Sin[7Pi*x]; ψ[x_]:=0; (*   initial conditions   *)

Plot[{ϕ[x],ψ[x]},{x,0,l}, PlotLabel -> "ϕ(x)"]

A[2]

A[3]

A[4]

A[5]

A[6]

A[7]

A[8]

(*   note that coefficients (Subscript[A, n])'s could be seen directly from ϕ(x) and that all (Subscript[B, n])'s are zero   *)

u[5,x,t]

Plot3D[Cos[5 π t] Sin[5 π x], {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 5])(x,t)"]   (* (Subscript[u, 5])(x,t) *)

u[7,x,t]

Plot3D[2 Cos[7 π t] Sin[7 π x], {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 7])(x,t)"]   (* (Subscript[u, 7])(x,t) *)

Plot3D[Cos[5 π t] Sin[5 π x]+2 Cos[7 π t] Sin[7 π x], {t,0,t0}, {x,0,l}, PlotLabel->"u(x,t)=(Subscript[(∑^∞), n = 1])(Subscript[u, n])(x,t) "]   (* solution u(x,t) *)

(*  Example 3. "The plucked string" *)

ϕ[x_]:= Piecewise[{{ 3x/2, x≤2/3}, { 3(1-x), x>2/3}}]; ψ[x_]:=0; (*   initial conditions   *)

Plot[{ϕ[x],ψ[x]},{x,0,l}, PlotLabel -> "ϕ(x)"]

A[1]

u[1,x,t]

Plot3D[((9 (Sqrt[3]) Cos[π t] Sin[π x])/(2 (π^2))), {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 1])(x,t)"]   (* (Subscript[u, 1])(x,t) *)

A[2]

u[2,x,t]

Plot3D[-((9 (Sqrt[3]) Cos[2 π t] Sin[2 π x])/(8 (π^2))), {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 2])(x,t)"]   (* (Subscript[u, 2])(x,t) *)

A[3]

A[4]

u[4,x,t]

Plot3D[((9 (Sqrt[3]) Cos[4 π t] Sin[4 π x])/(32 (π^2))), {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 4])(x,t)"]   (* (Subscript[u, 4])(x,t) *)

A[5]

u[5,x,t]

Plot3D[-((9 (Sqrt[3]) Cos[5 π t] Sin[5 π x])/(50 (π^2))), {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 5])(x,t)"]   (* (Subscript[u, 5])(x,t) *)

A[6]

A[7]

u[7,x,t]

Plot3D[((9 (Sqrt[3]) Cos[7 π t] Sin[7 π x])/(98 (π^2))), {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 7])(x,t)"]   (* (Subscript[u, 7])(x,t) *)

A[8]

u[8,x,t]

Plot3D[-((9 (Sqrt[3]) Cos[8 π t] Sin[8 π x])/(128 (π^2))), {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 8])(x,t)"]   (* (Subscript[u, 8])(x,t) *)

A[9]

A[10]

u[10,x,t]

Plot3D[((9 (Sqrt[3]) Cos[10 π t] Sin[10 π x])/(200 (π^2))), {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 10])(x,t)"]   (* (Subscript[u, 10])(x,t) *)

Plot3D[((9 (Sqrt[3]) Cos[π t] Sin[π x])/(2 (π^2)))-((9 (Sqrt[3]) Cos[2 π t] Sin[2 π x])/(8 (π^2)))+((9 (Sqrt[3]) Cos[4 π t] Sin[4 π x])/(32 (π^2)))-((9 (Sqrt[3]) Cos[5 π t] Sin[5 π x])/(50 (π^2)))+((9 (Sqrt[3]) Cos[7 π t] Sin[7 π x])/(98 (π^2)))-((9 (Sqrt[3]) Cos[8 π t] Sin[8 π x])/(128 (π^2)))+((9 (Sqrt[3]) Cos[10 π t] Sin[10 π x])/(200 (π^2))),{t,0,t0}, {x,0,l}, PlotLabel->"u(x,t)=(Subscript[(∑^∞), n = 1])(Subscript[u, n])(x,t) "]   (* solution u(x,t) *)

(*   Example 4. "Localized Plucking"   *)

p=2/3; a=p-0.2; b=p+0.1; h=0.5;

ϕ[x_]:= Piecewise[{{ 0, x≤a},{h*(x-a)/(p-a), a<x<p}, {h*(x-b)/(p-b), p<x<b}, { 0, x>b}}]; ψ[x_]:=0; (*   initial conditions   *)

Plot[{ϕ[x],ψ[x]},{x,0,l}, PlotLabel -> "ϕ(x)"]

A[1]

u[1,x,t]

Plot3D[0.134399Cos[π t] Sin[π x], {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 1])(x,t)"]   (* (Subscript[u, 1])(x,t) *)

A[2]

u[2,x,t]

Plot3D[-0.103467 Cos[2 π t] Sin[2 π x], {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 2])(x,t)"]   (* (Subscript[u, 2])(x,t) *)

A[3]

u[3,x,t]

Plot3D[-0.0375439 Cos[3 π t] Sin[3 π x], {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 3])(x,t)"]   (* (Subscript[u, 3])(x,t) *)

A[4]

u[4,x,t]

Plot3D[0.108307Cos[4 π t] Sin[4 π x], {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 4])(x,t)"]   (* (Subscript[u, 4])(x,t) *)

A[5]

u[5,x,t]

Plot3D[-0.0499331 Cos[5 π t] Sin[5 π x], {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 5])(x,t)"]   (* (Subscript[u, 5])(x,t) *)

A[6]

u[6,x,t]

Plot3D[-0.0350388Cos[6 π t] Sin[6 π x], {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 6])(x,t)"]   (* (Subscript[u, 6])(x,t) *)

A[7]

u[7,x,t]

Plot3D[0.0534346 Cos[7 π t] Sin[7 π x], {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 7])(x,t)"]   (* (Subscript[u, 7])(x,t) *)

A[8]

u[8,x,t]

Plot3D[-0.0211223 Cos[8 π t] Sin[8 π x], {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 8])(x,t)"]   (* (Subscript[u, 8])(x,t) *)

A[9]

u[9,x,t]

Plot3D[-0.00754167 Cos[9 π t] Sin[9 π x], {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 9])(x,t)"]   (* (Subscript[u, 9])(x,t) *)

A[10]

u[10,x,t]

Plot3D[0.0175493 Cos[10 π t] Sin[10 π x], {t,0,t0}, {x,0,l}, PlotLabel->"(Subscript[u, 10])(x,t)"]   (* (Subscript[u, 10])(x,t) *)

Plot3D[0.134399 Cos[π t] Sin[π x]-0.103467 Cos[2 π t] Sin[2 π x]-0.037544 Cos[3 π t] Sin[3 π x]+0.108307 Cos[4 π t] Sin[4 π x]-0.049933 Cos[5 π t] Sin[5 π x]-0.035039Cos[6 π t] Sin[6 π x]+0.053435 Cos[7 π t] Sin[7 π x]-0.021122 Cos[8 π t] Sin[8 π x]-0.007542 Cos[9 π t] Sin[9 π x]+0.017549 Cos[10 π t] Sin[10 π x],{t,0,t0}, {x,0,l}, PlotLabel->"u(x,t)=(Subscript[(∑^∞), n = 1])(Subscript[u, n])(x,t) "]   (* solution u(x,t) *)

(*  Hyperbolic Systems used in Secondary Oil Recovery*)
(*  Section 3 in the paper *)
(*  Consider the system D[(Subscript[s, w]),t]+D[(Subscript[f, w]),x]=0 and D[(Subscript[s, g]),t]+D[(Subscript[f, g]),x]=0  where (Subscript[s, w]) and (Subscript[s, g]) denote saturations of water and gas, (Subscript[s, o]):=1-(Subscript[s, w])-(Subscript[s, g]) is saturation of oil, and (Subscript[f, w]) and (Subscript[f, g]) denote the corresponding flux functions given as follows   *)
(*  Determine the eigenvalues of this system. Show that the system is hyperbolic since both eigenvalues are real. Show that there is only one umbilic point (i.e. point where the eigenvalues coincide)    *)
(Subscript[μ, w]):=0.5; (Subscript[μ, g]):=0.3; (Subscript[μ, o]):=1;

(Subscript[f, w])[sw_, sg_]:=((sw^2/(Subscript[μ, w]))/(sw^2/(Subscript[μ, w])+sg^2/(Subscript[μ, g])+(1-sw-sg)^2/(Subscript[μ, o])))

(Subscript[f, g])[sw_, sg_]:=((sg^2/(Subscript[μ, g]))/(sw^2/(Subscript[μ, w])+sg^2/(Subscript[μ, g])+(1-sw-sg)^2/(Subscript[μ, o])))

(* Jacobian Matrix*)

M={{D[(Subscript[f, w])[sw,sg], sw],D[(Subscript[f, w])[sw,sg], sg]},{D[(Subscript[f, g])[sw,sg], sw],D[(Subscript[f, g])[sw,sg], sg]}}

L=Eigenvalues[M]

(Subscript[λ, 1])[sw_,sg_]:=L[[1]]; (Subscript[λ, 2])[sw_,sg_]:=L[[2]];

Plot3D[L[[1]], {sw,0,1}, {sg,0,1}, PlotLabel-> "(Subscript[λ, 1])((Subscript[s, w]),(Subscript[s, g]))"]

Plot3D[L[[2]], {sw,0,1}, {sg,0,1}, PlotLabel-> "(Subscript[λ, 2])((Subscript[s, w]),(Subscript[s, g]))"]

Plot3D[L[[1]]-L[[2]], {sw,0,1}, {sg,0,1}, PlotLabel-> "(Subscript[λ, 1])-(Subscript[λ, 2])"]