Example 1.1: 2D compressible irrotational flow over a thin profile
The leading order correction of the velocity field due to the presence of the profile 
is governed by
is governed by
where 
is the Mach number.
is the Mach number.a) Determine the type of the system (3)
Let us re-write the system using matrices:
syms M
assume(M,'real'); assumeAlso(M>0);
A = [ 1-M^2, 0
0, 1 ];
B = [ 0, 1
-1, 0 ];
The characteristics can be computed by
syms dx dy
determinant = simplify( det( A'*dy - B'*dx ) )
Dx = simplify( solve( determinant, dx ) )
Therefore, 
would lead to two real characteristics and the system would thus be hyperbolic. On the other hand, no real characteristics can be found for 
and the flow would than be elliptic. The critical case 
is parabolic. This result is confirmed by experimental observation that supersonic flows contain shock waves, while subsonic flows are smooth.
would lead to two real characteristics and the system would thus be hyperbolic. On the other hand, no real characteristics can be found for and the flow would than be elliptic. The critical case is parabolic. This result is confirmed by experimental observation that supersonic flows contain shock waves, while subsonic flows are smooth.b) Convert (3) into a second-order PDE using velocity potential
Substituting the velocity potential 
, the equation (3b) is identically satisfied
, the equation (3b) is identically satisfied
such that only equation (3a) remains to be solved. After substitution:
We can classify second-order PDEs based on the coefficients multiplying the highest derivatives. The generic form of second order PDE reads
where in our case
a = [ 1-M^2, 0
0, 1 ];
b = [ 0
0 ];
The eigenvalues of a
eig(a)
have opposite sign if and (4) turns into a wave equation. On the other hand, the eigenvalues have the same sign if and (4) can be transformed to Poisson equation by stretching coordinates. The type of the problem is therefore conserved.
and (4) turns into a wave equation. On the other hand, the eigenvalues have the same sign if and (4) can be transformed to Poisson equation by stretching coordinates. The type of the problem is therefore conserved.Example 1.2 a): Find the characteristics if they exist
We know by now that characteristics exist for 
. The characteristics of a 2nd order PDE with 2 independent variables can be computed by solving
. The characteristics of a 2nd order PDE with 2 independent variables can be computed by solving
for 
. In our case
. In our case
where 
const. and 
const. defines the characteristics. We can then find the solution 
as a function of 
and 
, by converting the differentials
const. and const. defines the characteristics. We can then find the solution as a function of and , by converting the differentials
This brings the equation (4) to the normal form
with solution of the form
and 
can be determined by substituting the Ansatz (9) to boundary conditions. As further simplification, we use educated guess and say that shock waves should not travel upstream of the profile, therefore
where 
.
. Assume a parabolic profile
syms x y c epsilon
assume(x ,'real')
assume(y ,'real')
assume(c ,'real'); assumeAlso(c >0)
assume(epsilon,'real'); assumeAlso(epsilon>0)
p = 1-x^2;
y_p = epsilon*p;
p_x = diff(p,x)
The velocity on the profile surface must be tangential to the profile (no penetration boundary condition):
Subsitution of (10) to (11) gives
which after integration gives
where a remains an undetermined integration constant. Before and after the profile we have a symmetry plane at 
(no-penetration in case of irrotational flow):
(no-penetration in case of irrotational flow): 
which leads to
where b and d are integration constants. Due to the hyperbolicity we can assume that the incomming flow in front of the profile is undisturbed,
Moreover, we want a continuous velocity potential:
phi_0 = x;
xi = piecewise(y>=0,x-c*y,y<0,x+c*y)
phi_1 = piecewise(xi<=-1,0, xi>-1&xi<1,(xi^2-1)/c, xi>=1,0)
phi = phi_0 + epsilon*phi_1;
b)
Let's try to plug in some values
c = 2;
epsilon = 0.1;
u = diff(phi,x)
v = diff(phi,y)
and plot the solution
xmin=-2; xmax=2; ymin=-2; ymax=2; n=10;
figure(1); clf; hold on; axis([xmin xmax ymin ymax]);
fplot(subs(H),[-1 1]); fplot(-subs(H),[-1 1]);
fcontour(subs(phi),[xmin xmax ymin ymax]);
fcontour(subs(xi),'k--','LevelList',[-1 1]);
[x,y] = meshgrid(linspace(xmin,xmax,n),linspace(ymin,ymax,n));
quiver(x,y,subs(u),subs(v));
legend('upper surface','lower surface','velocity potential','characteristics','velocity field');
xlabel('x'); ylabel('y');
1.3: Time-dependent diffusion
a) Classification based on coefficients
The generic form of second order PDE reads
Let us select 
. Then the matrix 
and the vector 
for our case are given by
. Then the matrix and the vector for our case are given bysyms alpha
a = [ -alpha 0 0
0 -alpha 0
0 0 0 ];
b = [ 0
0
1 ];
The eigenvalues of
eig(a)
imply that the matrix is semi-definite, and the 
:
:rank([a b])
equals the number of independent variables. Thus we conclude that the PDE is parabolic.
b) Transformation to system of first-order equations
The original equation can be transformed into a system of first order equations
by introducing the auxiliary unknowns 
. In general we cannot say that the system is equivalent to the original equation, because the transformation can create spurious solutions. The solutions of (1) still solve (3), but not vice-versa. Therefore, the type of the first order system can also differ from the original equation. We will see that our example illustrates this problem. Let us write the system (3) in matrix form
. In general we cannot say that the system is equivalent to the original equation, because the transformation can create spurious solutions. The solutions of (1) still solve (3), but not vice-versa. Therefore, the type of the first order system can also differ from the original equation. We will see that our example illustrates this problem. Let us write the system (3) in matrix form
where 
and
andsyms alpha n_x n_y n_t
A = [ 1 0 0
0 0 0
0 0 0 ] ;
B = [ 0 -alpha 0
1 0 0
0 0 -1 ] ;
C = [ 0 0 -alpha
0 0 0
0 1 0 ] ;
We can analyze the system either with Fourier method, which analyzes whether solutions can have the form of a plane wave, or with the method of characteristics. Both methods lead to similar equations in the end. With the method of characteristics the normals to characteristic surfaces are given by
determinant = simplify( det( A'*n_t + B'*n_x + C'*n_y ) )
N_x = solve( determinant , n_x )
i.e. there is one real solution and two complex solutions. The real solution 
is typical for parabolic equations, since it represents the infinite speed of information spreading 
. The second component of this normal vector 
would be obtained by replacing the equation (3b) by
is typical for parabolic equations, since it represents the infinite speed of information spreading . The second component of this normal vector would be obtained by replacing the equation (3b) by
The two complex solutions 
most likely manifest the spurious solutions created by the non-equivalent transformation from (1) to (3). The system (3) is thus hybrid since it does not fall into any category. If we did not know the type of the original equation a priori, we would conclude that the original equation is not hyperbolic because it has for sure less than three real characteristics. One of the messages of this exercise is that second-order PDEs should be preferentially analysed by method a), unless it can be proved that the transformation to first-order system is equivalent. Two common cases where the transformation leads to equivalent first-order systems are considered in the next example. For more details refer to Wesseling (2001).
most likely manifest the spurious solutions created by the non-equivalent transformation from (1) to (3). The system (3) is thus hybrid since it does not fall into any category. If we did not know the type of the original equation a priori, we would conclude that the original equation is not hyperbolic because it has for sure less than three real characteristics. One of the messages of this exercise is that second-order PDEs should be preferentially analysed by method a), unless it can be proved that the transformation to first-order system is equivalent. Two common cases where the transformation leads to equivalent first-order systems are considered in the next example. For more details refer to Wesseling (2001).1.4 Steady 2D diffusion
Neglecting the variation of ϕ in time (corresponding to a steady/equilibrium state) leads to
The discriminant
implies the equation (7) is elliptic.
Using the same auxiliary unknowns as in 1.1b) we can eliminate ϕ from the transformed system to obtain the Cauchy-Riemann system
which is equivalent to the Laplace equation (7). Computation of characteristics
A = [ 1 0
0 -1 ];
B = [ 0 1
1 0 ];
determinant = simplify( det( A'*dy - B'*dx ) )
Dx = solve(determinant, dx)
does not reveal any real solutions, so the system is elliptic as well.
As we saw in 1.1), with more than two independent variables the transformation of both elliptic and parabolic equations creates spurious solutions and the resulting system is therefore no longer equivalent. Thus it is very often advantageous to analyze the type of an equaiton or a system of equations with several independent variables on sub-spaces of only two independent variables.
1.5 Korteweg-de Vries equation
We introduce the auxiliary unknowns 
such that we can form a system
such that we can form a system
Further we define a vector of unknowns 
and write the system (10) in matrix form as
and write the system (10) in matrix form as
where
A = [ 1 0 0
0 0 0
0 0 0 ];
B = [ 0 0 1
1 0 0
0 1 0 ];
The characteristics
determinant = simplify( det( A'*dx - B'*dt ) )
Dt = solve( determinant, dt )
are again typical for parabolic equations, representing infinite speed of information spreading.