Преобразование явного эйлера в неявный эйлер (через итерацию с фиксированной точкой)

Итак, у меня есть школьное задание, в котором мне нужно рассчитать положение группы автомобилей, следующих друг за другом на дороге (иначе говоря, движущихся по линии, поэтому, если Автомобиль 10 [автомобиль, который стоит первым] тормозит, тогда Автомобиль 9 тормозит). и когда Автомобиль 9 тормозит, Автомобиль 8 должен тормозить и т. д.).

Каждая машина следует правилу трех секунд (за исключением машины, которая стоит первой в очереди, она может ехать с любой скоростью, а все остальные машины в очереди соответствующим образом регулируют свою скорость). Каждая скорость автомобиля представлена ​​выражением:

введите здесь описание изображения

где «i» — это номер автомобиля, а «t» — момент времени (автомобиль с самым высоким индексом — это автомобиль, стоящий первым в очереди), а функция «f» представлена ​​следующим кодом:

% 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

Автомобиль, который находится впереди, просто имеет постоянную скорость 'g':

введите здесь описание изображения

В любом случае теперь, когда вы знаете контекст задачи, давайте перейдем к самой задаче. Это школьное задание состоит из нескольких шагов, и первым шагом является вычисление положения каждой машины во времени с использованием прямого/явного Эйлера, что я и сделал. Вот код для этого:

% 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

Однако проблема (где я застрял) заключается во втором шаге, где мне нужно теперь изменить свой код на первом шаге, чтобы использовать обратный/неявный эйлер через итерацию с фиксированной точкой. Я написал функцию, которая выполняет итерацию с фиксированной точкой, но помимо этого я действительно не знаю, что делать. Вот мой код для итерации с фиксированной точкой:

%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

Любая помощь приветствуется. Извините за длинный текст. Верхняя часть моего поста в основном представляет собой краткое изложение контекста задачи. Это не обязательно (и это не главное здесь), но я все равно добавил его, чтобы вы, ребята, знали, что происходит в моем коде.

Спасибо!


person Schytheron    schedule 08.02.2019    source источник


Ответы (1)


Проблема здесь в том, что вы написали свою итерацию с фиксированной точкой для скаляров, но у вас есть система дифференциальных уравнений. Если мы перепишем систему в векторной форме, она станет намного понятнее. Я добавил комментарий о том, как будут выглядеть явное и неявное уравнения таким образом, чтобы они действительно имели ту же стандартную форму y(t+1) = y(t) + h * f(y(t),t) (соответственно y(t+1) = y(t) + h * f(y(t+1),t+1)), которую вы можете найти в учебнике. Затем вы можете легко записать обновления:

% 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
person flawr    schedule 08.02.2019
comment
Большое спасибо! Я думаю, что понимаю общую суть этого, но я немного не знаком с некоторым синтаксисом MATLAB, который у вас есть (я довольно новичок в MATLAB), и я пытался его расшифровать, и я хочу прояснить некоторые вещи: 1. Я никогда не видел синтаксис [;] ([carFunc(posMat(c+1,t+1) - posMat(c,t+1)); g]) раньше, но, насколько я понимаю, вы создаете здесь матрицу, где результат слева от точки с запятой - это одна строка, а результат справа - другая строка... верно? продолжение ниже... - person Schytheron; 09.02.2019
comment
Да, квадратные скобки можно использовать для объединения матриц и векторов. Если вы используете запятые, вы получите горизонтальную конкатенацию, если вы используете точку с запятой, вы получите вертикальную конкатенацию. Ознакомьтесь с документами. - person flawr; 09.02.2019
comment
2. Можете ли вы объяснить, как работают последние два оператора «вывода» в «carFunc»? Что это должно значить? Как «выход» может принимать «маску» в качестве входного аргумента, если это не функция? - person Schytheron; 09.02.2019
comment
Доступ к матрицам MATLAB можно получить различными способами. Один из них называется логическое индексирование. . Это отлично подходит для доступа к записям на основе некоторых условий, и обычно это самый простой способ векторизации операторов if/else. Попробуйте, например, a = -5:5; a(a<0) = 0; a(a>2)=99;, чтобы понять. Он также работает в нескольких измерениях, как мы использовали его в carFunc. - person flawr; 09.02.2019
comment
3. Можно ли создавать циклы for без префикса for? Внутренний цикл for просто говорит c=1:M-1; без оператора for перед ним. Это просто трюк MATLAB, который работает и используется для сокращения кода? Извините, если надоедаю этими вопросами (и извините за дерьмовое форматирование). Я просто новичок в этом. - person Schytheron; 09.02.2019
comment
На самом деле это просто способ векторизации. Вы можете получить доступ ко всему диапазону элементов матрицы. Попробуйте например A = rand(5,5); A(1:2, 2:4) = 0;. В приведенной выше ссылке есть дополнительная информация об индексации MATLAB. - person flawr; 09.02.2019