468 pág.

# An Introduction to Programming and Numerical Methods in MATLAB

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