บทที่ 11
ความร้อนและสมการคลื่น
11.1
ส่วนนำ
การเอกสารเสริมการเรียนรู้นี้ออกแบบมาให้ผู้นสนใจได้คุ้นเคยกับบางวิธีการสำหรับแก้ปัญหาสมาการอนุพันธ์บางส่วนที่ขึ้นกับเวลา
โดยเฉพาะอย่างยิ่งเราจะนำวิธีการที่แตกต่างที่ชัดไปใช้ที่ได้นำเสนอไว้เป็นเบื้องต้นในการบรรยายเพื่อแก้สมการความร้อน
และสมการคลื่น ในเอกสารเสริมนี้จะนำผู้สนใจผ่านกระบวนการพัฒนาชุดของฟังก์ชันหนึ่ง
ซึ่งจะนำไปใช้ข้างหน้าและที่ผ่านมาของวิธีของ Euler สำหรับสมการความร้อนและวิธีกบกระโดดสำหรับสมการคลื่น เพื่อที่จะประหยัดเวลาในการพิมพ์ ผู้ใช้สามารถใช้วิธีคัดลอกรหัสคำสั่ง (code)
ที่อธิบายในเอกสารนี้มาอยู่ในไฟล์บนเครื่องซึ่งสามารถรันด้วยไซแลบ
โดยคิดให้ผู้ใช้ใช้ชื่อไฟล์ว่า tute_10.sce ที่ใส่รหัสคำสั่งทั้งหมดไว้
เพื่อเรียกสมาการความร้อนเข้าคือ
ut = cuxx 0 ≤ x ≤ 1,
t ≥ 0
u(0, x) = f(x)
u(t, 0) = α
u(t, 1) = β
และสมการคลื่นคือ
utt = cuxx 0 ≤ x ≤ 1,
t ≥ 0
u(0, x) = f(x)
ut(0, x) = g(x)
u(t, 0) = α
u(t, 1) = β.
11.2 สมการความร้อน
แรกสุดเราจะมุ่งเน้นไปที่สมการความร้อน วิธีง่ายที่สุดที่แตกต่างเด่นชัด อยู่บนฐานของวิธีของ Euler และกำหนดไว้ดังนี้
ที่ซึ่ง u j k j แทนคำตอบโดยประมาณ u(k△t, j△x) จากนี้เราทำให้ช่วง[0,1]เป็นดิสครีต เข้าไป/p n+1 points xj = j△x ที่ซึ่ง x = 1/n และช่วงเวลาเข้าไปในเวลาดิสครีต tk = k△t.
เพื่อนำไปใช้ประโยชน์ในตอนนี้ วิธี stepping เราต้องการเว็คเตอร์หนึ่งเพื่อบรรจุค่า u ไว้ที่ช่วงเวลาเฉพาะหนึ่ง ในความจริงเป็นประโยชน์ที่มีอีกเว็คเตอร์สำหรับเวลาที่ผ่านมา ขณะที่อินพุตให้กับฟังก์ชัน heat เราทำตัวบ่งชี้คือ
1. the number
of subintervals n in the x dimension,
2. the size of
the time step dt,
3. a function f
defining the initial conditions,
4. the left
boundary value alpha and
5. the right
boundary value beta.
ฟังก์ชันไซแลบของเราจะมีโครงสร้างพื้นฐานดังต่อไปนี้
//---------------------------------------------------------------
function [] =
heat(n,dt,f,alpha,beta)
//--------------------------
// Setup
parameters
//--------------------------
//--------------------------
// Setup
initial conditions
//--------------------------
//--------------------------
// Setup
Boundary conditions
//--------------------------
//--------------------------
// Plot the
initial conditions
//--------------------------
//--------------------------
// The main
time loop
//--------------------------
endfunction
//---------------------------------------------------------------
Cut and paste
this into your Scilab file heat.sce.
Now let us
expand the code. First, let us look at setting up parameters, as
this allows us
to define some variables for subsequent use. I chose the following
parameters
//--------------------------
// Setup
parameters
//--------------------------
c = 1.0;
dx = 1/n;
lambda =
dt/dx^2*c;
x = (0:dx:1)’;
t = 0;
Note that I am
using the comment lines to help you place the code correctly.
I needed c to
define the constant in the heat equation, dx is obvious, lambda
is a useful
combination of parameters that occurs in the updating of the solution
vector (by the
way, the dt is being passed in so that you can play with changing
it to test
stability). I have created an x vector so that I can plot our result x
versus u. A
variable for time t is useful for the main time loop. Cut and Paste
this code into
your file tute_10.sce in the initial part of the function.
The initial
condition will be given by the value of a function f at the x values.
It turns out
that we also need a second u vector u1 which we create at this point.
We simply use
//--------------------------
// Setup
initial conditions
//--------------------------
u0 = f(x);
u1 = zeros(u0);
We will need to
define a function f. Lets do something simple like a sin curve.
I used
//---------------------------------------------------------------
function
u=fsin(x)
k=1
u =
sin(%pi*k*x)
endfunction
//---------------------------------------------------------------
Computational
Science, Scilab Tutorials 75
Put this in
another file util.sce file. This file will hold extra functions which
can be used
independently of the heat equation function.
The boundary
conditions can be simply set to α at the left hand end of u0
and β at
the right hand end.
//--------------------------
// Setup
Boundary conditions
//--------------------------
u0(1) = alpha;
u0($) = beta;
In most
simulations it is useful to plot out results. This is particularly useful
when coding up
a method as it often points to errors in the code. So lets plot
the initial
condition. I suggest something like
//--------------------------
// Plot the
initial conditions
//--------------------------
xselect()
myplot(x,u0)
The xselect()
will bring the graphics window to the front and the myplot
function is
just a function to plot x versus u vectors. Copy this code to the
correct
position within your file.
Now we also
have to define our plotting routine. I used
//---------------------------------------------------------------
function
[]=myplot(x,u)
wait = 50000;
xbasc();
plot2d(x,u,rect=[0,-1,1,1]);
xpause(wait)
endfunction
//---------------------------------------------------------------
The main part
of this function is just the standard plot2d command. I have
added the
xpause to make the plots look nicer when we have our time evolution
working. The
command leaves the plot up for at least the specified number of
microseconds.
Put a copy of this function in your utility file util.sce.
Lets try our
code. Save your file and in the Scilab command window load
in the
functions and run the heat function.
exec heat.sce
exec util.sce
n=20;
heat(n,0.5/n^2,fsin,0,-1)
Computational
Science, Scilab Tutorials 76
The right
boundary condition doesn’t match up with the initial conditions.
We could run
our method with this incompatibility, but let’s fix up the boundary
conditions.
With new boundary conditions try
n=20;
heat(n,0.5/n^2,fsin,0,0)
Nothing
happens.
Now we can go
back the code and add in the time evolution part. We need a
loop, we can
use a for loop but I like using a while loop based on time. I.e. I
suggest using
//--------------------------
// The main
time loop
//--------------------------
while t<1
t=t+dt
//--------------------------------
// Update the
solution vector
//--------------------------------
myplot(x,u0)
end
This loop is a
little incorrect, in that it may overshoot the final time by dt.
Not to worry.
The main thing now is how to update the solution vector. The
most obvious
way is to use a for loop to implement the finite difference update
uk+1
j = uk
j +
c_t
_x2 uk
j+1 − 2
uk
j + uk
j−1_
The uk
j values will
be stored in u0 and the new values uk+1
j in the vector
u1.
Once the update
is completed we will over write u0 with the new u1. The code
is
//--------------------------------
// Update the
solution vector
//--------------------------------
for j=2:n
u1(j) =
(1-2*lambda)*u0(j) ...
+
lambda*u0(j+1) + lambda*u0(j-1);
end
u0 = u1;
Computational
Science, Scilab Tutorials 77
It is always
important to get the range of j correct. We only need to work
with the
interior points. There are actually n+1 entries in u0 and u1 and so the
range j=2:n
corresponds to interior points.
An alternative
method to the for loop is to use the array notation in Scilab
//--------------------------------
// Update the
solution vector
//--------------------------------
u1(2:$-1) =
(1-2*lambda)*u0(2:$-1) ...
+
lambda*u0(3:$) + lambda*u0(1:$-2);
u0 = u1;
This second
method will be more efficient in Scilab though the for method
is closer to
what you would use in C or Fortran.
Perhaps use the
first method initially and change over to the second if the
evolution seems
slow. As we will be plotting the solution for every iteration the
overall speed
of the method will be slow anyway. The array method becomes
more important
for higher dimensional problems.
Lets see if our
method works. Add in the time loop code, load in the code
with exec and
run the heat command again. You might like to increase the
number of grid
points to say 100.
To stop the
evolution you can use the Control menu with the pause, resume
and abort
commands.
Lets try
another initial condition. Go back to your code and add the following
function to
your utility file.
//---------------------------------------------------------------
function
u=fline(x)
u = 0.5-x;
endfunction
//---------------------------------------------------------------
This is a
solution of the Heat equation with the boundary conditions α =
0.5
and β = −0.5.
The exact solution to the heat equation with this set of boundary
conditions is
the straight line between these two points. Let’s check how our
solver works.
exec heat.sce;
exec util.sce
n=100;heat(n,0.5/n^2,fline,0.5,-0.5)
Why is the
solution changing? Go back to your code and find the bug and
fix it so that
the boundary condition is set properly.
Computational
Science, Scilab Tutorials 78
11.2.1
Stability
Our analysis of
the method tells us that we need _t < _x2
2c for
stability. Let’s
see what
happens if we go over this constraint. Try
n=100;
heat(n,0.52/n^2,fsin,0.0,0.0)
After some time
you should see instabilities arising.
11.3 Wave
Equation
Let’s quickly
produce an implementation of the leap frog method for the wave
equation. We
need to implement the algorithm
uk+1
j = 2 uk
j − uk−1
j + c
_t2
_x2 uk
j+1 − uk
j + uk
j−1_
We will use
three vectors u2 to hold uk+1
j , u1 to hold
uk
j and u0 to
hold uk−1
j .
The structure
is similar to that for the heat equation.
//---------------------------------------------------------------
function [] =
wave(n,dt,f,g,alpha,beta)
//--------------------------
// Setup
parameters
//--------------------------
c = 1.0;
dx = 1/n;
lambda =
dt^2/dx^2*c;
x = (0:dx:1)’;
t=0;
//--------------------------
// Setup
initial and
// boundary
conditions
//--------------------------
u0 = f(x);
u0(1) = alpha;
u0($) = beta;
u1 = zeros(u0);
if (g==1) then
u1(2:$-1) =
u0(1:$-2);
elseif (g==-1)
then
Computational
Science, Scilab Tutorials 79
u1(2:$-1) =
u0(3:$);
else
u1(2:$-1) =
0.5*(u0(1:$-2)+u0(3:$));
end
u1(1) = alpha;
u1($) = beta;
u2 = zeros(u0);
u2(1) = alpha;
u2($) = beta;
//--------------------------
// Plot the
initial conditions
//--------------------------
xselect()
myplot(x,u0)
myplot(x,u1)
//--------------------------
// The main
time loop
//--------------------------
while t<20
t=t+dt;
//--------------------------------
// Update the
solution vector
//--------------------------------
for j=2:n
u2(j) =
(2-2*lambda)*u1(j) - u0(j) ...
+
lambda*u1(j+1) + lambda*u1(j-1);
end
u0=u1;
u1=u2;
myplot(x,u2)
end
endfunction
//---------------------------------------------------------------
11.3.1
Experimentation
We can use the
initial conditions for u that we used for the heat equation.
Add in the code
in the previous subsection to a wave.sce file. Load in the
new functions
using exec. Now play around with the code. For instance try
Computational
Science, Scilab Tutorials 80
n=100;
wave(n,1/n,fsin,0,0,0)
Hopefully you
will see an “oscillating string”.
Lets try a
different initial condition, A bump in the middle of the string.
//---------------------------------------------------------------
function u =
fbump(x)
c = 0.5; w =
0.125;
cl = c-w;
cr = c+w;
u =
bool2s(x>cl & x<=cr)
endfunction
//---------------------------------------------------------------
See how we can
use comparisons to create interesting functions.
Now try out
this initial condition
n=100; wave(n,1/n,fbump,0,0,0)
It looks cool
as the disturbance reflects off the boundaries. It is amazing, we
are actually
getting the exact solution.
You can get the
bump to move either to the right or left using the g argument.
Try
n=100;
wave(n,1/n,fbump,1,0,0)
and
n=100;
wave(n,1/n,fbump,-1,0,0)
11.3.2
Stability
Try testing the
stability of the method for slightly larger timesteps. For even
slightly larger
timesteps than 1/n you should see growing modes in the solution,
especially if
you use the bump as an initial condition. The effect is not seen so
quickly if you
use smoother initial conditions.
11.4 Exercises
1. You might
like to play with a solver for the heat equation which uses the
implicit Euler
method. Here is the code
Computational
Science, Scilab Tutorials 81
//---------------------------------------------------------------
function [] =
iheat(n,dt,f,alpha,beta)
//--------------------------
// Setup
parameters
//--------------------------
c = 1.0;
dx = 1/n;
lambda =
dt/dx^2*c;
x = (0:dx:1)’;
t=0;
//--------------------------
// Call f to
setup initial
// conditions.
u1 set to same
// size as u0
//--------------------------
u0 = f(x);
//--------------------------
// Setup the
matrix A as a
// sparse
matrix
//--------------------------
A = (1+2*lambda)*diag(sparse(ones(n-1,1)))
...
-lambda*diag(sparse(ones(n-2,1)),-1)
...
-lambda*diag(sparse(ones(n-2,1)),1);
//--------------------------
// Setup
Boundary conditions
//--------------------------
u0(1) = alpha
u0($) = beta
//--------------------------
// Plot the
initial conditions
//--------------------------
xselect()
myplot(x,u0)
//--------------------------
// The main
time loop
Computational
Science, Scilab Tutorials 82
//--------------------------
while t<20
t = t+dt
u0(2:$-1) = A\(u0(2:$-1))
myplot(x,u0)
end
endfunction
//---------------------------------------------------------------
Try
n=100;
iheat(n,1/n,fsin,0,0)
It evolves so
quickly you might want to change the wait time in the myplot
function.
Try
n=100;
iheat(n,1/n,fbump,0,0)
2. Use the code
for heat and iheat to produce a program that solves the heat
equation using
the Crank Nicolson method. Try a smooth initial condi-
tion like fsin
and a rough initial condition like fbump. Something strange
happens with
the rough initial condition.
You might like
to experiment with Crank Nicolson using fsin, but with
progressively
high frequencies. Can you make a conjecture about the rela-
tive efficiency
of Crank Nicolson to damp out errors at different frequencies.
How does this
explain the observed behaviour of Crank Nicolson method
on the fbump initial data.
