matlabfixed-point-iteration

Converting explicit euler to implicit euler (via fixed-point iteration)


So I have a school task where I need to calculate the position of a bunch of cars following each other on a road (AKA driving in a line, so if Car 10 [the car that is first in line] brakes, then Car 9 brakes and when Car 9 brakes, Car 8 has to brake etc.).

Each car follows the "three second rule" (except the car that is first in line, he can drive as fast as he wants and all other cars in line adjust their velocity accordingly). Each cars velocity is represented by the expression:

enter image description here

where 'i' is the index of the car and 't' is a point in time (the car with the highest index is the car that is first in line) and the function 'f' is represented by this code:

% Input x: Distance to the car in front (of you. NOT the car that is at the
% very front of the line).
% Output: Car velocity
function output = carFunc(x)
    if (x >= 75)
        output = 25;
    else
        output = x/3;
    end
end

The car that is at the very front just has a constant velocity 'g':

enter image description here

Anyway, now that you know the context of the task lets actually move on to the task itself. This school task consists of multiple steps and the first step is to calculate the position of each car over time using forward/explicit Euler, which I have done. Here is the code for that:

% Euler's Method
% Initial conditions and setup
g = 25; % Velocity of car M (the car that is first in line, ahead of all others) [m/s]
h = 0.01;  % step size
M = 10; % Number of cars
x = 0:h:40;  % the range of x (time)
n = numel(x);  % the number of timesteps
posMat = zeros(M,n); % Matrix that holds positions of each car (car = row, time = column)
for i=1:M
    posMat(i,1) = 10*i;  % the initial positions for the cars
end
for t=1:(n-1)
    % Calculate position of all cars EXCEPT car M (using the first
    % equation)
    for c=1:(M-1)
    f = carFunc(posMat(c+1,t) - posMat(c,t)); % Velocity of car 'c'
    posMat(c,t+1) = posMat(c,t) + h * f; % Explicit Euler
    end
    % Calculate positon of last car M (first car in line) here:
    % x_M^(n+1) = x_M^n + h*g
    posMat(M,t+1) = posMat(M,t) + h*g;
end

The problem however (where I got stuck) is the second step where I need to now modify my code in the first step to use backwards/implicit euler via fixed-point iteration. I have written a function that does fixed-point iteration but beyond that I don't really know what to do. Here is my code for fixed-point iteration:

%Fixed-Point Iteration
%Computes approximate solution of g(x)=x
%Input: function handle g, starting guess x0,
% number of iteration steps k
%Output: Approximate solution

function out=fpi(g, x0, k)
    x = zeros(1, k+1);
    x(1)=x0;
        for i=1:k
            x(i+1)=g(x(i));
        end
    out=x(k+1);
end

Any help is appreciated. Sorry for the long text. The top part of my post is mainly just a "short" summary of the context behind the task. It isn't absolutely necessary (and it isn't the focus here) but I added it anyway so you guys know what is going on in my code.

Thanks!


Solution

  • The problem here is that you wrote your fixed point iteration for scalars but you have a system of differential equations. If we rewrite the system in vectorial form it becomes much clearer. I added a comment on how what the explicit and the implicit equation would look like that way such that they really have the same standard form y(t+1) = y(t) + h * f(y(t),t) (resp. y(t+1) = y(t) + h * f(y(t+1),t+1)) that you can find in a texbook. Then you can easily write down the updates:

    % Euler's Method
    % Initial conditions and setup
    g = 25; % Velocity of car M (the car that is first in line, ahead of all others) [m/s]
    h = 0.01;  % step size
    M = 10; % Number of cars
    x = 0:h:40;  % the range of x (time)
    n = numel(x);  % the number of timesteps
    posMat = zeros(M,n); % Matrix that holds positions of each car (car = row, time = column)
    k = 20; % number of fixed point iterations
    explicit = true;
    for i=1:M
        posMat(i,1) = 10*i;  % the initial positions for the cars
    end
    
    for t=1:n-1 
        % Calculate position of all cars EXCEPT car M (using the first
        % equation)
        c=1:M-1;
        if explicit
            f = [carFunc(posMat(c+1,t) - posMat(c,t)); g]; % Velocity of car 'c'
            posMat(:,t+1) = posMat(:,t) + h * f; % Explicit Euler
        else 
            %explicit euler:
            %posMat(:,t+1) = posMat(:,t) + h * [carFunc(posMat(c+1,t) - posMat(c,t)); g];
            %implicit euler: (needs to be solved)
            %posMat(:,t+1) = posMat(:,t) + h * [carFunc(posMat(c+1,t+1) - posMat(c,t+1)); g];
    
            %fixed point iteration:
            posMat(:,t+1) = posMat(:,t); % initialization
            for m=1:k
                posMat(:,t+1) = posMat(:,t) + h * [carFunc(posMat(c+1,t+1) - posMat(c,t+1)); g]; 
            end
        end
    end
    plot(posMat');
    
    
    % Input x: Distance to the car in front (of you. NOT the car that is at the
    % very front of the line).
    % Output: Car velocity
    function output = carFunc(x)
        mask = x >= 75;
        output = zeros(size(x));
        output(mask) = 25;
        output(~mask) = x(~mask)/3;
    end