ME 17 - Homework #5 Solving Partial Differential Equations Poisson's Equation
ME 17 - Homework #5 Solving Partial Differential Equations Poisson's Equation
ME 17 - Homework #5 Solving Partial Differential Equations Poisson's Equation
Richard Behiel
ME 17 Gibou
June 6
th
, 2014
ME 17 Homework #5
Solving Partial Differential Equations
Poissons Equation
2
The Assignment
Consider a domain , a given function
Suppose that we sample the solution at n points.
(a) Write down the general form of the linear system associated with the standard
discretization for
.
(b) Write an algorithm that builds the matrix and the right-hand-side of the linear system
given as well as the boundary conditions.
(c) Solve this system using the \ operator in Matlab
(d) Compare your results with the exact solution for
and . Also try the biggest problem you can solve on your computer, i.e.
take as big as you can without crashing your machine.
(e) Write a report explaining in details every step above and upload it on Gauchospace
along with your Matlab code. The report can include some figures comparing the
numerical solution to the exact solution.
2. (Extra Credit) Consider the advection equation
where and Show that, for this equation, the numerical method
is always unstable. Here and are the steps in time and space, respectively.
3
Introduction to Poissons Equation
Poissons equation is an elliptic partial differential equation which is commonly used to model diffusion
in a wide range of fields, such as electrostatics, mechanical engineering, and theoretical physics. Named
after French scientist and mathematician Simon Denis Poisson, the equation is expressed as follows:
where and are functions on a manifold. In one-dimensional Cartesian coordinates, where and
are functions of , the above equation is identical to
If were given by a known analytic function, then could be solved for in the above equation by
integrating twice and applying boundary conditions to the solution. However, in practice is often a
set of measured points without any simple analytic representation, so discrete numerical methods must
be used in order to solve for . The domain of and will be discretized with constant spatial
increment such that the domain is evenly divided into grid points, where is some arbitrarily large
integer. The discrete method used to numerically solve Poissons equation relies on the finite difference
approximation of the second derivative,
where a subscript on a function represents the function being evaluated at that particular grid point,
and is some integer in [ ]. Applying Poissons equation to the left hand side of the above equation,
Multiplying by
Writing the LHS as a vector product,
[
][
This equation relates and in the vicinity of the
, where
[
,
[
and
=
[
The rows of A, , and
, so it can
be plotted against the discretized domain vector as an approximation of the function ; no further math
is required. An example plot of the function specified in the assignment, with , is shown in
Figure 1. Entries of are connected by straight lines. The error for such a low value of is clearly
visible, but is negligible for sufficiently large .
Figure 1: Plot of and for
5
The Matlab Code
In this section is the Matlab function that I have written and uploaded with this report, which takes as
input and outputs a plot of the discrete solution in blue, as well as the exact function in red. It does
this by constructing matrix A (variable A in the code) and vector
have been initialized, the function solves for with the backslash command and plots the results.
Every step of the code is commented adequately.
%%%ME 17 Homework #5%%%
%Solves the discrete Poisson's equation with n spatial intervals
function Poisson(n)
%%%Parameters%%%
x1 = 0; x2 = 1; %Define boundary conditions
u_x1 = 0; u_x2 = 0; %Define initial conditions
x = linspace(x1,x2,n); %x values, domain vector
dx = x(2) - x(1); %x increment, delta x
dx2 = dx*dx; %Simplifies calculation by squaring before
f = @(x) -4*(pi^2)*sin(2*pi*x); %Definition of function f
%%%Matrix A and b vector%%%
A = sparse(n,n); %Initialize matrix A
b = zeros(n,1); %Initialize b vector
A(1,1) = 1; b(1) = u_x1; %Apply B.C. and initialize A(1,1)
A(n,n) = 1; b(n) = u_x2; %Apply B.C. and initialize A(n,n)
for i = 2:n-1 %Loop through b and diagonal of A
b(i) = f(x(i))*dx2; %Set value i of b to f_i given by equation
A(i,i) = -2; %Diagonal values are all -2
A(i,i-1) = 1; %Values on either side of diagonal are 1, as
A(i,i+1) = 1; %shown in the matrix A
end %End of for loop
%%%Plotting the results
u = A \ b; %Solve for discrete u vector
x_ = linspace(x1,x2,100); %Set domain for exact u
plot(x_,sin(2*pi*x_),'r') %Plot exact u onto the same plot
hold on %Keep plot open
plot(x,u,'b') %Plot discrete u vector
xlabel('x') %Label x axis
ylabel('u(x)') %Label u axis
title(['Discrete Poisson Equation, n = ',num2str(n)]); %Give plot a title
hold off %Plot can no longer be added to
end %End of function
This function uses the sparse command when initializing the matrix A, which decreases the amount of
memory required to run the function. Matrix A contains
Attempting to solve this equation with the equation numerical method shown below will always lead to an
unstable solution.
To see why, solve the equation for
)
Then, dividing the equation by
)
Now, rearranging the equation,
At this point, the equation can be interpreted in such a way that we can tie the behavior of some arbitrary
over space to its behavior over time. The leftmost term of the equation is the scaling ratio of
as it goes
from time step to . If this ratio is greater than 1, then
)
Since and are positive, this condition can be reduced to
Moving terms around,
Multiplying by the denominator,
10
The above equation is a spatial condition which is tied to the behavior of
grows from one grid point in space to the next,
also grows from the current time step to the next. This, in
turn, will make
even greater than it was previously, and the cycle repeats. A similar argument can show
that the same concept applies to