[Math] Implement a program in Matlab for LU decomposition with pivoting

computer sciencemath-softwareMATLABnumerical methods

I need to write a program to solve matrix equations Ax=b where A is an nxn matrix, and b is a vector with n entries using LU decomposition. Unfortunately I'm not allowed to use any prewritten codes in Matlab. I am having problems with the first part of my code where i decompose the matrix in to an upper and lower matrix.

I have written the following code, and cannot work out why it is giving me all zeros on the diagonal in my lower matrix. Any improvements would be greatly appreciated.

function[L R]=LR2(A)

%Decomposition of Matrix AA: A = L R

z=size(A,1);

L=zeros(z,z);

R=zeros(z,z);

for i=1:z

% Finding L

for k=1:i-1

L(i,k)=A(i,k);

for j=1:k-1

L(i,k)= L(i,k)-L(i,j)*R(j,k);

end

L(i,k) = L(i,k)/R(k,k);

end

% Finding R

for k=i:z

R(i,k) = A(i,k);

for j=1:i-1

R(i,k)= R(i,k)-L(i,j)*R(j,k);

end

end

end

R

L

end

I know that i could simply assign all diagonal components of L to be 1, but would like to understand what the problem is with my program!

I am also wondering how to change this program to include pivoting. I understand I need to say that if a diagonal element is equal to zero something needs to be changed. How would I go about this?
Thanks in advance!

Best Answer

For backward and forward elimination I used

% Now use a vector y to solve 'Ly=b' 
 for j=1:z

    for k=1:j-1
       b(j)=b(j)-L(j,k)*b(k);
    end;
    b(j) =b(j)/L(j,j); 
 end;

% Now we use this y to solve Rx = y

x = zeros( z, 1 );

for i=z:-1:1    
 x(i) = ( b(i) - R(i, :)*x )/R(i, i); 
end

I put that into your code, and it works ;) For helping you with pivot strategies it would be helpful to know what strategie you want/have to use as there are various versions. One simple version would be to just swap rows such that the diagonal element $a_{ii} \neq 0$ for all $i = 1, \dots, z$.

Related Solutions

[Math] A problem for Matlab solver

You are looking for a shooting method.

In other words, you have the following problem description: $$\mathbf{y}'(t) = f(\mathbf{y},t),$$ $$y_1(0) = 0,$$ $$y_2(\tau) = 0.$$

This is a boundary value problem, and you need to transform it into an initial value problem.

To do so, instead of solving the equations above, you are really solving $$\mathbf{y}'(t) = f(\mathbf{y},t),$$ $$y_1(0) = 0,$$ $$y_2(0) = b.$$

for some estimate of $b$. Your problem is to find $b$ such that the original condition is met.

One technique is known as shooting. You make an initial guess for $b$, compute the error at the right boundary for $y_2(\tau)$, and adjust $b$ based on your error. This is essentially a root-finding problem, and now you can use root-finding techniques to find the "root" $b$.

Essentially, the solution to an ODE initial value problem can be written as $y(t;y_0)$, or the solution $y(t)$ parameterized by the initial conditions $y_0$. This is because, broadly speaking, solutions to ODEs are not unique until you specify initial conditions. So the actual solution trajectory depends greatly on your choice of initial conditions $y_0$.

In your case, you have a second-order system, so you can compute $\mathbf{y}(t;\mathbf{y}_0)$, where $\mathbf{y}_0 = \begin{pmatrix} 0 \\ b\end{pmatrix}$. Since the trajectory that leads you to $y_2(\tau) = 0$ is unique, you can iterate over different values of $b$ to find the right result, and indeed, formulate the problem as solving $F(b) = 0$, where $F$ is your "function" that computes the solution to your ODE. In other words, treat the MATLAB code that you've written above as a function of the parameter $b$ that returns the value $y_2(\tau)$. Then, compute the root.

There are a number of ways to do this; the simplest method, although by far not the fastest, is through bisection. Newton's method may also work, but that is more difficult to set up. Broyden's method (the secant method) will also work.

EDIT: Below is a solution. If this is homework, you must attribute this code properly. This code is presented to show you how shooting works, not to solve the problem for you

The following code can be pasted in a single file called "odetest.m". This is not the highest-quality code; I banged it together to illustrate a point, and I tried to design it in such a way so as to be mathematical in how it was structured. Certainly it can be optimized. That's not the point.

First, we define an RK4 ODE solver. Then, we define a function $F(b)$ such that $F$ returns the value of $y_2(t=1)$. This value depends entirely on the parameter $b$; as such, we wish to find $b$ such that $F(b) = 0$.

We then use Broyden's method, which is a secant method, to compute the estimate for $b$ from an initial guess of both the value and the derivative. Finally, we call this, and then plot some things.

The resulting plot is here: Broyden's Method Shooting

Broyden's method is not the only way to solve this problem; however, it converges quite quickly and as you can see, it is easy to implement.

function odetest
    [i,lambda,lh] = broyden(@F,0,3);

    figure;
    hold on
    colors = 'brgky';
    legend_string = {};
    for j=1:length(lh)
        F(lh(j),1,colors(mod(j-1,length(colors))+1));
        legend_string{j} = ['j  = ' num2str(j) ', b = ' num2str(lh(j))];
    end
    legend(gca,legend_string);

    plot([0 1],[0 0],'k'); % Plot the x-axis

    set(gcf,'PaperPositionMode','auto')
    print('-f1','-r360','-dpng','broydenshooting.png');

    answer = lambda
end

function x=F(b,ploton,color)
    % Extra arguments for plotting mode; when the function is called as a
    % function, just use a single argument
    ode = @(t,y)[-y(1)+y(2);-1+y(2)];
    y0 = [0; b];
    T = [0 1];
    h = 0.05;
    t = T(1):h:T(2);

    y = RK4(ode,.05,t,y0);
    if nargin > 1
        plot(t,y(2,:),color,'LineWidth',2)
    end
    x = y(2,end);
end

function y=RK4(f,h,t,y0)
    n = length(t);
    y = zeros(length(y0),n);
    y(:,1) = y0;
    for i=1:n-1
        k1 = f(t(i),y(:,i));
        k2 = f(t(i)+h/2, y(:,i)+h/2*k1);
        k3 = f(t(i)+h/2, y(:,i)+h/2*k2);
        k4 = f(t(i)+h, y(:,i)+h*k3);
        y(:,i+1) = y(:,i)+h/6*(k1+2*k2+2*k3+k4);
    end
end

function [i,lambda,lh] = broyden(func,IX,IJ)
% lambda = broyden(func,IX,IJ)
% IX -- initial guess for X
% IJ -- initial guess for J = dF/dX
%
% i -- number of steps needed to converge
% lambda -- solution
% lh -- lambda history, the computed lambda at each step
%
% Attempt to compute the solution to a system of nonlinear equations
% without explicitly requiring the computation of the derivatives of the
% function.
lh = [];
lambda = IX;
B = IJ;
    for i=1:30
        lh = [lh; lambda];
        f = func(lambda);
        s = B\(-f);
        x = s + lambda;
        y = func(x)-f;
        B = B + ((y-B*s)*s')/(s'*s);
        if norm(x-lambda)/norm(lambda) < 1e-10
            lambda = x;
            break
        end
        lambda = x;
    end
end

[Math] Strictly column diagonally dominant matrices and Gaussian elimination with partial pivoting

It seems to me that the following computation proves the result.

Let ${{a'}}_{j , k}$ be the coefficients after the first step i.e. for $j > 1$

$${{a'}}_{j , k} = {a}_{j , k}-{a}_{j , 1} \frac{{a}_{1 , k}}{{a}_{1 , 1}}$$

Let us compute

$$\renewcommand{\arraystretch}{2} \begin{array}{rcl}\displaystyle \sum\limits _{\substack{j = 2\\j \neq k}}^{n} \left|{{a'}}_{j , k}\right|&=&\displaystyle \sum\limits _{\substack{j = 2\\j \neq k}}^{n} \left|{a}_{j , k}-{a}_{j , 1} \frac{{a}_{1 , k}}{{a}_{1 , 1}}\right|\\ & \leqslant &\displaystyle \sum\limits _{\substack{j = 2\\j \neq k}}^{n} \left|{a}_{j , k}\right|+\sum\limits _{\substack{j = 2\\j \neq k}}^{n} \left|{a}_{j , 1}\right| \frac{\left|{a}_{1 , k}\right|}{\left|{a}_{1 , 1}\right|}\\ & < &\displaystyle \left|{a}_{k , k}\right|-\left|{a}_{1 , k}\right|+\left|{a}_{1 , k}\right| \sum\limits _{\substack{j = 2\\j \neq k}}^{n} \frac{\left|{a}_{j , 1}\right|}{\left|{a}_{1 , 1}\right|}\\ & \leqslant &\displaystyle \left|{a}_{k , k}\right|-\left|{a}_{1 , k}\right|+\left|{a}_{1 , k}\right| \left(1-\frac{\left|{a}_{k , 1}\right|}{\left|{a}_{1 , 1}\right|}\right)\\ & \leqslant &\displaystyle \left|{a}_{k , k}\right|-\left|{a}_{1 , k}\right| \frac{\left|{a}_{k , 1}\right|}{\left|{a}_{1 , 1}\right|}\\ & \leqslant &\displaystyle \left|{a}_{k , k}-{a}_{1 , k} \frac{{a}_{k , 1}}{{a}_{1 , 1}}\right|\\ & \leqslant &\displaystyle \left|{{a'}}_{k , k}\right| \end{array}$$

Best Answer

Related Solutions

[Math] A problem for Matlab solver

[Math] Strictly column diagonally dominant matrices and Gaussian elimination with partial pivoting

Related Question