NotebookFileName[]

Preliminaries

This is just to show how we make a string from several pieces. Below I use Epilog-> to place to line segments on the graph.

Plot[1+Sin[x],{x,0,Pi},PlotRange->{{-1+0,Pi+1},{-2,2}},Epilog->{{Line[{{-1,0},{0,1}}],Line[{{Pi,1},{Pi+1,0}}]}},AspectRatio->(1/5),ImageSize->600]

Below I show how to “animate” the above graph using Cos[t] as a time function. The point is that I have to adjust the line segments as well.

Manipulate[Plot[(1+Sin[x])Cos[t],{x,0,Pi},PlotRange->{{-1+0,Pi+1},{-2,2}},Epilog->{{Line[{{-1,0},{0,1*Cos[t]}}],Line[{{Pi,1*Cos[t]},{Pi+1,0}}]}},AspectRatio->(1/5),ImageSize->600],{t,0,2 Pi}]

Statement of the Problem

Below is the PDE + BCs + ICs that we want to solve here.

PDE: (∂_tt)u[x,t]=(∂_xx)u[x,t], 0≤x≤π, t≥0, 

BCs:  u[0,t]=0, u[π,t]+(∂_x)u[π,t]=0,

ICs:u[x,0]=ff1[x], (∂_t)u[x,0]=0.

For demonstration purposes, we will look for the positive eigenvalues only and their corresponding eigenfunctions:

The general solution can be obtained using Mathematica function DSolve

Solving the Boundary-Eigenvalue Problem

Finding the eigenvalues

DSolve[-D[A[x],{x,2}]==(μ^2) A[x],A[x],x]

Below I show how to substitute the special value x=0 in the general solution. Mathematica code /. is called the substitution rule.

((C[1]) Cos[x μ]+(C[2]) Sin[x μ])/.{x->0}

The second boundary condition is more complicated. I substitute the length of the string Pi.

((C[1]) Cos[x μ]+(C[2]) Sin[x μ]+D[(C[1]) Cos[x μ]+(C[2]) Sin[x μ],{x,1}])/.{x->Pi}

I use the Collect function to group terms with C[1] and C[2]. (Here Mathematica displays C[1] and C[2] in a nicer for in Output.)

Collect[((C[1]) Cos[x μ]+(C[2]) Sin[x μ]+D[(C[1]) Cos[x μ]+(C[2]) Sin[x μ],{x,1}])/.{x->Pi},{C[1],C[2]}]

Thus, the system that we have to solve for C[1] and C[2] is:

(C[1])=0
(C[2]) (μ Cos[π μ]+Sin[π μ])+(C[1]) (Cos[π μ]-μ Sin[π μ])=0

Of course we can solve this system “by hand”, but I do it slow. The determinant of the system is: (Below, I write the matrix and its determinant how Mathematica does it: A matrix is in Mathematica a list of its rows.)

Det[
{ (* the matrix starts here *)
{1,0}, (* the first row, coefficients only *)
{(Cos[π μ]-μ Sin[π μ]),(μ Cos[π μ]+Sin[π μ])} (* the second row, coefficients only *)
}(* the matrix starts here *)
]

Our system has a nontrivial solution for C[1] and C[2] if and only if

μ Cos[π μ]+Sin[π μ]=0

Let us see on a graph where the solutions are:

Plot[μ Cos[π μ]+Sin[π μ],{μ,0,10}]

The command solve cannot find the exact solutions of the equation:

μ Cos[π μ]+Sin[π μ]=0

Try

Solve[μ Cos[π μ]+Sin[π μ]==0,μ]

The command that will find approximate solutions is FindRoot

FindRoot[μ Cos[π μ]+Sin[π μ]==0,{μ,0.5}]

The logic of this command is that we give an initial guess for μ and Mathematica finds a very accurate approximate value.

Below is the FindRoot as a function of the initial guess:

(μ/.FindRoot[μ Cos[π μ]+Sin[π μ]==0,{μ,#}])&[0.5]

Verify that the value found is accurate solution:

μ Cos[π μ]+Sin[π μ]/.{μ-> 0.7876372941648639`}

Quite close to 0.

Now try other solutions of

μ Cos[π μ]+Sin[π μ]=0

Plot[μ Cos[π μ]+Sin[π μ],{μ,0,10}]

One close to 2

(μ/.FindRoot[μ Cos[π μ]+Sin[π μ]==0,{μ,#}])&[2]

Now try more, say ten, and name the table of values:

mymus=Table[(μ/.FindRoot[μ Cos[π μ]+Sin[π μ]==0,{μ,#}])&[s],{s,0.5,10.5,1}]

Verify on the graph if these are good solutions. The values that we calculated are at the red points.

Map[Point[{#,0}]&,mymus]

Plot[μ Cos[π μ]+Sin[π μ],{μ,0,11},Epilog->{{RGBColor[1,0,0],Map[Point[{#,0}]&,mymus]}}]

The values seem quite accurate. But we can try more of them

Clear[mymus];
mymus=Table[(μ/.FindRoot[μ Cos[π μ]+Sin[π μ]==0,{μ,#}])&[s],{s,0.5,20.5,1}]

Verify on the graph if these are good solutions. The values that we calculated are at the red points.

Plot[μ Cos[π μ]+Sin[π μ],{μ,0,21},Epilog->{{RGBColor[1,0,0],Map[Point[{#,0}]&,mymus]}}]

Let us verify the values numerically as well. Here we use Map[] function

(# Cos[π #]+ Sin[π #])&

is so called Pure Function, instead of the variable we put # and finish the Pure Function with &

Map[(# Cos[π #]+ Sin[π #])&,mymus]

The above is quite impressive. We can make it even more impressive by Chop[]

Chop[Map[(# Cos[π #]+ Sin[π #])&,mymus]]

Finding the eigenfunctions

Thus, we have 21 values for μ-s which are in a list which we call mymus; the squares of μ-s  are the eigenvalues. Now write the corresponding eigenfunctions. For the μ-s that we found we need a nontrivial solution of the system

(C[1])=0
(C[2]) (μ Cos[π μ]+Sin[π μ])+(C[1]) (Cos[π μ]-μ Sin[π μ])=0

A nontrivial solution is C[1]=0 and C[2]=1. Thus the corresponding eigenfunctions are: (Go back to the general solution and plug in C[1]=0 and C[2]=1 and use mus in the list mymus:

((C[1]) Cos[x μ]+(C[2]) Sin[x μ])/.{C[1]->0,C[2]->1}

Clear[myeigfuns];myeigfuns[x_]=Sin[x μ]/.{μ->mymus}

The eigenfunctions are in the list called myeigfuns[x]. The Mathematica code for selecting particular elements of a list is as follows, the second element

myeigfuns[x][[2]]

These are our eigenfunctions. I made a claim that they are orthogonal to each other. Let us verify that, just say the first two

myeigfuns[x][[1]]

myeigfuns[x][[2]]

Show that they are orthogonal:

Integrate[myeigfuns[x][[1]] myeigfuns[x][[2]],{x,0,Pi}]

Now try the entire table of all integrals. That will be 21x21 matrix which should have only diagonal elements nonzero: (The command below is slow. Using NIntegrate, numerical integration would be faster. But, NIntegrate might require special options to be successful. You can find the options in the bigger notebook that I posted on the class website.)

21 21

Table[Chop[Integrate[myeigfuns[x][[j]] myeigfuns[x][[k]],{x,0,Pi}]],{j,1,21},{k,1,21}]

Constructing natural modes of vibrations

Next is to do time part. That is similar. Time part has two functions Cos and Sin. Let us ignore Sin, which corresponds to the initial velocity. The time Cos parts are

Cos[t μ]/.{μ->mymus}

Clear[mytimeCs];mytimeCs[t_]=Cos[t μ]/.{μ->mymus}

mytimeCs[t][[1]]

Now we have found the natural modes of vibrations of this string: The first natural mode of vibration is

myeigfuns[x][[1]]mytimeCs[t][[1]]

See the value at at t=0

myeigfuns[x][[1]](mytimeCs[0][[1]])

Plot[myeigfuns[x][[1]]mytimeCs[0][[1]],{x,0,Pi},PlotRange->{{0,Pi+1},{-1,1}},Epilog->{{Line[{{Pi,myeigfuns[Pi][[1]]mytimeCs[0][[1]]},{Pi+1,0}}]}},AspectRatio->(1/5),ImageSize->600]

Let us animate it:

Manipulate[Plot[myeigfuns[x][[1]]mytimeCs[t][[1]],{x,0,Pi},PlotRange->{{0,Pi+1},{-1,1}},Epilog->{{Line[{{Pi,myeigfuns[Pi][[1]]mytimeCs[t][[1]]},{Pi+1,0}}]}},AspectRatio->(1/5),ImageSize->600],{t,0,((2 Pi)/(mymus[[1]]))}]

Manipulate[Plot[myeigfuns[x][[2]]mytimeCs[t][[2]],{x,0,Pi},PlotRange->{{0,Pi+1},{-1.1,1.1}},Epilog->{{Line[{{Pi,myeigfuns[Pi][[2]]mytimeCs[t][[2]]},{Pi+1,0}}]}},AspectRatio->(1/5),ImageSize->600],{t,0,((2 Pi)/(mymus[[2]]))}]

Manipulate[Plot[myeigfuns[x][[k]]mytimeCs[t][[k]],{x,0,Pi},PlotRange->{{0,Pi+1},{-1.1,1.1}},Epilog->{{Line[{{Pi,myeigfuns[Pi][[k]]mytimeCs[t][[k]]},{Pi+1,0}}]}},AspectRatio->(1/5),ImageSize->600],{t,0,((2 Pi)/(mymus[[k]])),Appearance->{"Open"}},{k,Range[21],Setter,Appearance->{"Open"}}]

Solving the Initial Value Problem

Now introduce a specific initial condition. It must satisfy the boundary conditions. Let us try with a quadratic function (x^2)+ b x and choose b to satisfy the second boundary condition:

(((x^2)+ b x)/.{x->Pi})+(D[((x^2)+ b x),x])/.{x->Pi}

Solve[((((x^2)+ b x)/.{x->Pi})+(D[((x^2)+ b x),x])/.{x->Pi})==0,b]

Try the following initial condition

((x^2)+ b x)/.{b->((-2 π-(π^2))/(1+π))}

Clear[ff];ff[x_]=(1/4)((((-2 π-(π^2)) x)/(1+π))+(x^2))

Plot[ff[x],{x,0,Pi},PlotRange->{{0,Pi+1},{-1,1}},Epilog->{{Line[{{Pi,ff[Pi]},{Pi+1,0}}]}},AspectRatio->(1/5),ImageSize->600]

myeigfuns[x][[4]]

Integrate[ff[x]myeigfuns[x][[1]],{x,0,Pi}]/Integrate[myeigfuns[x][[1]]myeigfuns[x][[1]],{x,0,Pi}]

Clear[coeffs];coeffs=Table[Integrate[ff[x]myeigfuns[x][[k]],{x,0,Pi}]/Integrate[myeigfuns[x][[k]]myeigfuns[x][[k]],{x,0,Pi}],{k,1,21}]

Clear[appff];appff[x_]=Sum[coeffs[[k]] myeigfuns[x][[k]],{k,1,21}]

Plot[{ff[x],appff[x]},{x,0,Pi},PlotRange->{{0,Pi},{-1,1}},AspectRatio->(1/5),ImageSize->600]

An approximation for the solution is:

Clear[appuu];appuu[x_,t_]=Sum[coeffs[[k]] myeigfuns[x][[k]] mytimeCs[t][[k]],{k,1,21}]

Test it with specific values:

appuu[1,1]

Animate

Manipulate[Plot[appuu[x,t],{x,0,Pi},PlotRange->{{0,Pi+1},{-1,1}},Epilog->{{Line[{{Pi,appuu[Pi,t]},{Pi+1,0}}]}},AspectRatio->(1/5),ImageSize->600],{t,0,8,Appearance->{"Open"}}]

Solving the Initial Value Problem with a more Complicated f(x)

In this section I introduce a more complicated function f(x) for which symbolic Integrate[] can not calculate the coefficients. We use NIntegrate[] with specific options.

I just multiply the preceding initial displacement with a high power of the Sin[] function:

Clear[ff1];ff1[x_]=(1/4)((((-2 π-(π^2)) x)/(1+π))+(x^2))((Sin[x])^6)

Satisfies the boundary conditions:

ff1[x]/.{x->0}

(ff1[x]/.{x->Pi})-D[ff1[x],{x,1}]/.{x->Pi}

See the function

Plot[ff1[x],{x,0,Pi},PlotRange->{{0,Pi+1},{-1,1}},Epilog->{{Line[{{Pi,ff[Pi]},{Pi+1,0}}]}},AspectRatio->(1/5),ImageSize->600]

myeigfuns[x][[4]]

See how Integrate[] is slower than NIntegrate[]

tI=Timing[Integrate[ff1[x]myeigfuns[x][[1]],{x,0,Pi}]/Integrate[myeigfuns[x][[1]]myeigfuns[x][[1]],{x,0,Pi}]]

Pay attention how we need special options in NIntegrate[]

tNI=Timing[NIntegrate[ff1[x]myeigfuns[x][[1]],{x,0,Pi},MaxRecursion->200,AccuracyGoal->10,PrecisionGoal->10]/NIntegrate[myeigfuns[x][[1]]myeigfuns[x][[1]],{x,0,Pi},MaxRecursion->200,AccuracyGoal->10,PrecisionGoal->10]]

NIntegrate is

((tI[[1]])/(tNI[[1]]))

times faster.

Below I calculate the coefficients for the initial coefficient ff1[x]

Clear[coeff1s];coeff1s=Table[NIntegrate[ff1[x]myeigfuns[x][[k]],{x,0,Pi},MaxRecursion->200,AccuracyGoal->10,PrecisionGoal->10]/NIntegrate[myeigfuns[x][[k]]myeigfuns[x][[k]],{x,0,Pi},MaxRecursion->200,AccuracyGoal->10,PrecisionGoal->10],{k,1,21}]

Here is our approximation of the initial condition ff1[x]

Clear[appff1];appff1[x_]=Sum[coeff1s[[k]] myeigfuns[x][[k]],{k,1,21}]

How good is it?

Plot[{ff1[x],appff1[x]},{x,0,Pi},PlotStyle->{{RGBColor[0,0,0.7],Thickness[0.01]},{RGBColor[1,0,0],Thickness[0.004]}},PlotRange->{{0,Pi},{-1,1}},AspectRatio->(1/5),ImageSize->600]

Quite good approximation!

An approximation for the solution of PDE+BCs+ICs:

Clear[appuu1];appuu1[x_,t_]=Sum[coeff1s[[k]] myeigfuns[x][[k]] mytimeCs[t][[k]],{k,1,21}]

Test it with specific values:

appuu1[1,1]

Animate

Manipulate[Plot[appuu1[x,t],{x,0,Pi},PlotRange->{{0,Pi+1},{-1,1}},Epilog->{{Line[{{Pi,appuu1[Pi,t]},{Pi+1,0}}]}},AspectRatio->(1/5),ImageSize->600],{t,0,12.8,Appearance->{"Open"}}]