Reija Ristilä, 43177L, 29.1. 1999

Problem 1.1

The classical numbering scheme leads to block-tridiagonal matrix A (with n=4):

        [ B   I   0   0 ]                [-4   1   0   0 ]
   A =  [ I   B   I   0 ]    where  B =  [ 1  -4   1   0 ]
        [ 0   I   B   I ]                [ 0   1  -4   1 ]
        [ 0   0   I   B ]                [ 0   0   1  -4 ]

According to Kreyszig, in matrix A there are blocks I, not -I, as is written in the lecture notes. The point of the exercise is to assemble matrix A using Matlab from blocks or diagonals.

Using blocks

n = 4;                                 % can be varied
N = n^2;
B = 4*diag(ones(1,n))-diag(ones(1,n-1),1)-diag(ones(1,n-1),-1);
O = zeros(n,n);
pI= eye(n,n);

first = [ B pI zeros(n,(N-2*n)) ];     % first row  
middle = [ pI B pI zeros(n,(N-3*n))];  % building block for rows after the 1st

A = first;
for i = 0:(n-2)
  A = [A; zeros(n,i*n) middle(:,1:(N-i*n))]
end
spy(A)

Using diagonals

function A = lapm(n)
% LAPM(n) assembles matrix A for 2- dim problem u''= 0.
% n is the block size, n^2 is the number of unknown nodes

N = n^2;
d1 = ~~mod(1:N,n)';                  % the (1 1 ..1 0 1 1 ..) diagonal
d2 = ones(N,1);
A = spdiags([d2 d1],[-n -1],N,N);          % lower triangle matrix
A = A + A' + spdiags(-4*ones(N,1),0,N,N);

First only a lower triangular matrix is created to avoid problems with the function spdiags.

Problem 1.2

Solving of the steady-state temperature leads to linear equation Ax=-B, where A is the Laplacian and in B there are the boundary values. Reshaping the boundary conditions into a vector:

left=[20,40,60,80];top=[50,20,10,0];bot=zeros(1,4);right=left;
m=5;n=m-1;
B=[left' zeros(n,(n-2)) right']+[bot; zeros((n-2),n); top];
B=flipud(B)';B=B(:);

Creating matrix A and solving the equation Ax=-B, then reshaping the solution

A=lapm(n);
Tvek= A\-B;
T=reshape(Tvek,4,4)';   % reshape works columnwise, thus you need to transpose

Adding boundary conditions:

T=[fliplr(left)' T fliplr(right)'];  
NW=0.5*(left(n)+top(1));NE=0.5*(top(n)+right(n));  
SE=0.5*(bot(n)+right(1));SW=0.5*(bot(1)+left(1));
T=[NW top NE; T; SW bot SE];
surfc(T), colorbar

graph The result can be checked by calculating the differences:

for i=2:n; for j=2:n ;
  tt(i,j)= 4*T(i,j)-T(i-1,j)-T(i+1,j)-T(i,j-1)-T(i,j+1);
end; end

tt =
   1.0e-13 *
         0         0         0         0
         0   -0.0711    0.1066         0
         0   -0.1421   -0.1776    0.1421
         0   -0.1066   -0.0711   -0.0355

which is very close to zero.

Problem 1.3

function U = lapsolve1(left,top,right,bot)
% U = lapsolve1(left,top,right,bot);
% Solves Laplace equation on a square. Arguments are
% the boundary values, each of length n. 
% U is the solution matrix with the boundaries (n+2)x(n+2).

n = max(size(left));
B = [left' zeros(n,(n-2)) right']+[bot; zeros((n-2),n); top];
B = flipud(B)';B = B(:); 

A = lapm(n);
Tvek = A\-B;
T = reshape(Tvek,n,n)';
  
T = [fliplr(left)' T fliplr(right)'];  
NW = 0.5*(left(n)+top(1)); NE = 0.5*(top(n)+right(n));  
SE = 0.5*(bot(n)+right(1)); SW = 0.5*(bot(1)+left(1));
T = [NW top NE; T; SW bot SE];

Testing the function with a finer grid:

x = [ 1 2 3 4]; x2 = 1:0.5:4;
p = polyfit(x,left,2); left2 = polyval(p,x2);
p = polyfit(x,top,2); top2 = polyval(p,x2);
p = polyfit(x,right,2); right2 = polyval(p,x2);
p = polyfit(x,bot,2); bot2 = polyval(p,x2);

U = lapsolve1(left2,top2,right2,bot2);

Checking by calculating the differences, which are again small:

tt =
   1.0e-13 *
         0         0         0         0         0         0         0
         0   -0.1421         0    0.0711    0.1066   -0.0355   -0.0711
         0   -0.0711    0.0711   -0.1066   -0.0355   -0.0711         0
         0    0.0711   -0.0711   -0.2842    0.1776   -0.2842   -0.1421
         0   -0.0711    0.1066   -0.0355   -0.1066         0   -0.2842
         0    0.1066    0.1066   -0.2487         0   -0.1421         0
         0   -0.0355   -0.0355    0.1421    0.0355   -0.1421    0.0355

graph

Problem 1.4

x = 1:0.5:12;
top = 100*ones(size(x));
left = zeros(size(x));
right = left; bot = left;

U = lapsolve1(left,top,right,bot);
surfc(flipud(U)),colorbar

graph

Problem 1.5

A) Equation for a damped pendelum:

v'' + cv' + k sin(v) = 0

To reduce it to a first-order system, it can be written as:

[ y1 ] = v    =>   / y1' = y2
[ y2 ] = v'        \ y2' = -k sin(y1) -c y2

The critical points of the system are found at y' = 0, which leads to (y1 = +/- n*pi, y2 = 0).

The system can be linearized at (0,0) by approximating sin(x) with x. Hence

 y' = [  0   1 ] y
      [ -k  -c ]

For other points, a change of variables must be done to be able to linearize the system. This leads to change in the sign of element A(2,1), that is, -k => k.

B) Eigenvalues for the first system are:

b1 = 1/2 ( -c + sqrt(c^2 -4k) ), b2 = 1/2 ( -c - sqrt(c^2 -4k) )

and for the second system:

bb1 = 1/2 ( -c + sqrt(c^2 +4k) ), bb2 = 1/2 ( -c - sqrt(c^2 +4k) )

For determining the types of the critical points, for the first system it applies:

p = b1 + b2 = -c,   q = b1 * b2 = k,   det = p^2 - 4q = c^2 -4k

and for the second system:

p = -c, q = -k,  det = c^2 +4k

The type of the points depends on the values of the constants c and k such that a critical point is

node if q > 0 and det >= 0,
saddle if q < 0,
center if p = 0 and q > 0
spiral if p <> 0 and det < 0.

C) Odefile for solving the system of differential equations:

function dy = ex15(t,y,flag,k,c)
dy =[
y(2);
-c*y(2)-k*sin(y(1))];

graphs

Reija Ristilä

Last modified: 1999-02-03