An Introduction to Programming and Numerical Methods in MATLAB
468 pág.

An Introduction to Programming and Numerical Methods in MATLAB


DisciplinaMatlab498 materiais2.121 seguidores
Pré-visualização50 páginas
dsolve(\u2019Dy = t\u2c62+y\u2c62\u2019,\u2019y(0)=0\u2019,\u2019t\u2019);
N = 20;
tv = linspace(eps,1,N);
for i = 1:N
f(i) = double(subs(y,t,tv(i)));
end
Here we have used the commands
dsolve Which solves the di\ufb00erential equation de\ufb01ned as its \ufb01rst argument,
subject to the initial condition y(0)=0 and \ufb01nally we use the third argu-
ment to tell MATLAB that t is the independent variable.
8.3 Banded Matrices 259
subs This substitutes values from the array tv in the solution of the di\ufb00erential
equation.
double This returns a double precision value (that is it converts a symbolic
value to a double precision value).
We shall now discuss a class of methods which can be derived using Taylor
series.
8.3 Banded Matrices
In many problems matrices will not only be sparse but they will also have very
well-de\ufb01ned structure. Let us consider a tri-diagonal matrix, so that the only
non-zero entries are in the super- and sub-diagonal. We consider the set of
equations
a1y2 + b1y1 = r1,
a2y3 + b2y2 + c2y1 = r2,
... =
...
an\u22121yn + bn\u22121yn\u22121 + cn\u22121yn\u22122 = rn\u22121,
bnyn + cnyn\u22121 = rn.
Notice that we can solve the \ufb01rst equation to give y1 in terms of y2, so that
y1 =
r1
b1
\u2212 a1
b1
y2
which can then be substituted into the next equation to eliminate y1. This gives
a2y3 + b2y2 + c2
(
r1
b1
\u2212 a1
b1
y2
)
= r2,
a2y3 +
(
b2 \u2212 c2 a1
b1
)
y2 = r2 \u2212 c2 r1
b1
.
This amounts to a rede\ufb01nition of b2 and r2. This equation can then be used to
express y2 in terms of y3 and then substitute it into the next equation, until
we reach the \ufb01nal equation at which point we have an equation solely in yn
which is easily solved. The system at this point has only got a super-diagonal
(we have eliminated the sub-diagonal) and hence the values of all the preceding
y\u2019s can be found by back substitution. This is written in code as:
260 8. Solving Di\ufb00erential Equations
\ufffd
\ufffd
\ufffd
\ufffd
% Set up system
x = 0.0:pi/6:pi;
N = length(x);
h = x(2)-x(1);
a = 1/h\u2c62*ones(size(x));
b = -2/h\u2c62*ones(size(x));
c = 1/h\u2c62*ones(size(x));
r = x.*sin(x);
a(1) = 0; b(1) = 1; r(1) = 0;
c(N) = 0; b(N) = 1; r(N) = 0;
% Forward sweep
for j = 2:N
b(j) = b(j) - c(j)*a(j-1)/b(j-1);
r(j) = r(j) - c(j)*r(j-1)/b(j-1);
end
% Final equation
y(N) = r(N)/b(N);
for j = (N-1):-1:1
y(j) = r(j)/b(j)-a(j)*y(j+1)/b(j);
end
xf = 0:pi/50:pi;
sol = 2*(1-cos(xf))-xf.*sin(xf)-4*xf/pi;
clf
plot(x,y,\u2019.\u2019,\u2019MarkerSize\u2019,24)
hold on
plot(xf,sol,\u2019b\u2019)
xlabel(\u2019x\u2019,\u2019FontSize\u2019,15)
ylabel(\u2019f\u2019,\u2019FontSize\u2019,15)
text(0.5,-0.2,\u2019Finite difference solution\u2019,\u2019FontSize\u2019,12)
text(0.7,-0.3,\u2019shown as blobs\u2019,\u2019FontSize\u2019,12)
This is actually solving the following problem
d2y
dx2
= x sinx y(0) = y(\u3c0) = 0.
We can solve this by hand by integrating twice and applying the boundary
conditions to obtain
y(x) = \u22122 cosx \u2212 x sinx \u2212 4x
\u3c0
+ 2.
8.3 Banded Matrices 261
As you can see even using 7 points determines the solution quite well
This method is called a Thomas algorithm.
Example 8.5 In order to solve the system
d2y
dx2
+ 4
dy
dx
\u2212 y = cosx
subject to the boundary conditions y(0) = 0 and y(\u3c0) = 1, we use the code
0 0.5 1 1.5 2 2.5 3 3.5
\u22121.8
\u22121.6
\u22121.4
\u22121.2
\u22121
\u22120.8
\u22120.6
\u22120.4
\u22120.2
0
x
f
Finite difference solution
shown as blobs
262 8. Solving Di\ufb00erential Equations
\ufffd
\ufffd
\ufffd
\ufffd
% Set up system
x = 0.0:pi/60:pi;
N = length(x);
h = x(2)-x(1);
a = (1/h\u2c62+4/(2*h))*ones(size(x));
b = (-2/h\u2c62-1)*ones(size(x));
c = (1/h\u2c62-4/(2*h))*ones(size(x));
r = cos(x);
a(1) = 0; b(1) = 1; r(1) = 0;
c(N) = 0; b(N) = 1; r(N) = 1;
% Forward sweep
for j = 2:N
b(j) = b(j) - c(j)*a(j-1)/b(j-1);
r(j) = r(j) - c(j)*r(j-1)/b(j-1);
end
% Final equation
y(N) = r(N)/b(N);
for j = (N-1):-1:1
y(j) = r(j)/b(j)-a(j)*y(j+1)/b(j);
end
clf
plot(x,y)
xlabel(\u2019x\u2019,\u2019FontSize\u2019,15)
ylabel(\u2019y\u2019,\u2019FontSize\u2019,15)
which gives
8.4 Runge\u2013Kutta Methods 263
0 0.5 1 1.5 2 2.5 3 3.5
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x
y
We have used the discretisation
yn+1 \u2212 2yn + yn\u22121
h2
+ 4
yn+1 \u2212 yn\u22121
2h
\u2212 yn = cosxn.
This can be extended to penta-diagonal systems, which can again be done
using a Thomas algorithm but this becomes more and more complicated. How-
ever, we can also exploit the sparse command. We can also solve periodic
problems which give rise to sparse problems, but with elements in the top right
and bottom left corners in addition to the diagonals.
8.4 Runge\u2013Kutta Methods
We shall not derive these methods in fullness but shall give a \ufb02avour of how
this may be done. In order to solve the system
dy
dt
= f(t, y),
we use Taylor series and as such we write the scheme as
yn+1 = yn + ak1 + bk2 (8.3)
where
k1 = \u2206t f(tn, yn),
k2 = \u2206t f(tn + \u3b1\u2206t, yn + \u3b2k1).
264 8. Solving Di\ufb00erential Equations
In order to derive an equation of this form we start with the Taylor series for
yn+1 = y(t + \u2206t)
yn+1 = yn + \u2206tf(tn, yn) +
\u2206t2
2
f \u2032(tn, yn) + · · · ,
where the prime denotes a derivative with respect to t. Using the chain rule we
can write
df
dt
=
\u2202f
\u2202t
+
\u2202f
\u2202y
dy
dt
,
but using the original equation (y\u2d9 = f)
df
dt
=
\u2202f
\u2202t
+
\u2202f
\u2202y
f,
so that
yn+1 = yn + \u2206tf(tn, yn) +
\u2206t2
2
(
\u2202f
\u2202t
+
\u2202f
\u2202y
f
)
· · · .
Let us now substitute the forms for k1 and k2 into Equation (8.3), so that
yn+1 = yn + a\u2206tf(tn, yn) + b\u2206tf (tn + \u3b1\u2206t, yn + \u3b2\u2206tf(tn, yn)) .
We now need to expand the last term in this equation, so that
f (tn + \u3b1\u2206t, yn + \u3b2\u2206tf(tn, yn)) \u2248 f(tn, yn) + \u3b1\u2206t \u2202f
\u2202t
\u2223\u2223\u2223\u2223
tn,yn
+ \u3b2\u2206tf(tn, yn)
\u2202f
\u2202y
\u2223\u2223\u2223\u2223
tn,yn
.
By substitution and comparison with the other Taylor series we \ufb01nd that
in \u2206tf a + b = 1,
in \u2206t2
\u2202f
\u2202t
\u3b1b =
1
2
,
in \u2206t2f
\u2202f
\u2202y
\u3b2b =
1
2
.
We have only three equations but with four unknowns, so we are free to choose
one value. For instance a = 1/2 implies b = 1/2 with in turn gives \u3b1 = \u3b2 = 1.
The scheme is then
yn+1 = yn +
\u2206t
2
(f(tn, yn) + f (tn + \u2206t, yn + \u2206tf(tn, yn))) .
We can of course choose other values of a (or any of the other variables). In
order to code this algorithm we use:
8.4 Runge\u2013Kutta Methods 265
\ufffd
\ufffd
\ufffd
\ufffd
%
% runge_kutta.m
%
t = 0.0:0.1:5.0;
del_t = 0.1;
a = 0.5; b = 1-a;
alpha = 0.5/b; beta = 0.5/b;
n = length(t);
y = zeros(size(t));
y(1) = 0;
for ii = 1:n-1
time = t(ii);
k_1 = del_t*func(t(ii),y(ii));
k_2 = del_t*func(t(ii)+alpha*del_t,y(ii)+beta*k_1);
y(ii+1) = y(ii) + a * k_1 + b * k_2;
end
exact = exp(t-t.\u2c62/2);
using
\ufffd
\ufffd
\ufffd
function [value] = func1(t,y)
value = y*(1-t);
Using this number of points there is no visible di\ufb00erence between the numerical
and the exact solution. We can investigate the e\ufb00ect of altering the value of a.
This method can be extended to higher orders: for instance the fourth-order
scheme is given by
yn+1 = yn +
1
6
(k1 + 2k2 + 2k3 + k4)
where
k1 = \u2206tf(tn, yn)
k2 = \u2206tf
(
tn +
\u2206t
2
, yn +
k1
2
)
k3 = \u2206tf
(
tn +
\u2206t
2
, yn +
k2
2
)
k4 = \u2206tf (tn + \u2206t, yn + k3) .
MATLAB o\ufb00ers two commands ode23 and ode45 which perform these in-
tegrations. It uses a more advanced numerical scheme which iterates to correct
the solution as it progresses. The above code could be replaced by
266 8. Solving Di\ufb00erential Equations
\ufffd
\ufffd
\ufffd
\ufffd
y0 = 1; tspan = [0 5];
[tt,yy] = ode45(\u2019func1\u2019,tspan,y0);
This uses an adaptive step which is based on the local gradients. Interestingly
this chooses to use 57 points to attain the default accuracy (which is de\ufb01ned
in terms of relative and absolute errors, see the manual pages help ode45 for
details).
8.5 Higher-Order Systems
8.5.1 Second-Order Systems