My numerical solution does not align with my exact solution

I am trying to minimized this problem

where

and

Exact solutions are

My matlab code is below(note: had some of the code from the book we are using);

    global oldu oldt_lambda u_new x1 x2 t_lambda tfwd_interp 
   
%parameters
    delta = .001;
    N =1000;                  
    tin = linspace(0, 1, N+1)'; 
    tfwd_interp = tin;
    t_lambda = tin;
    
    % Initialization
    u_new = zeros(N+1, 1);      
    x1 = zeros(N+1, 1);         
    x2 = zeros(N+1, 1);         
    lambda1 = zeros(N+1, 1);   
    lambda2 = zeros(N+1, 1);    
    t_bwd = 0;

function dydt = state_equations(t, y)
    global u_new t_lambda 

    u_interp = interp1(t_lambda, u_new, t, 'spline'); 
 
    dx1dt = y(2,:);    
    dx2dt = u_interp;  
  
    dydt = [dx1dt; dx2dt];
end

function dldt = costate_equations(t,l)
     
    dlambda1_dt = 0;  
    dlambda2_dt = -1*l(1,:)-1;   
    
    dldt = [dlambda1_dt; dlambda2_dt];
end

    
    test = -1;
    j = 1;

    while (test < 0 && j < 100)
        oldtfwd = tfwd_interp;
        oldt_lambda = t_lambda;
        oldu = u_new;
        oldx1 = x1;
        oldx2 = x2;
        oldlambda1 = lambda1;
        oldlambda2 = lambda2;
        
        % Forward 
        tspan = linspace(0, 1, N+1);
        y0 = [0; 0];  
        [t_fwd, yret] = ode45(@(t, y) state_equations(t, y), tspan, y0);
        x1 = yret(:,1);
        x2 = yret(:,2);
        tfwd_interp = t_fwd;
        
        % Backward
        lambda_tspan = linspace(1, 0, N+1);
        lambda0 = [0; 0];  
        [t_bwd, Lambda] = ode45(@(t, l) costate_equations(t, l), lambda_tspan, lambda0);
        
        lambda1 = flip(Lambda(:, 1));
        lambda2 = flip(Lambda(:, 2));
        %t_lambda = flip(t_bwd);
        
        % Control update
        t_lambda=t_bwd;
        u1 = -0.5*lambda2;  
        u_new = 0.5 * (u1 + oldu);  
        
        % Convergence check
        if j > 1
            temp1 = delta * sum(abs(u_new)) - sum(abs(oldu - u_new));
            temp2 = delta * sum(abs(x1)) - sum(abs(oldx1 - x1));
            temp3 = delta * sum(abs(x2)) - sum(abs(oldx2 - x2));
            temp4 = delta * sum(abs(lambda1)) - sum(abs(oldlambda1 - lambda1));
            temp5 = delta * sum(abs(lambda2)) - sum(abs(oldlambda2 - lambda2));
           % if (temp1 >= 0 && temp2 >= 0 && temp3 >= 0 && temp4 >= 0 && temp5 >= 0)
                %test = 0;
            %end
        end
        j = j + 1;
    end

    % exact 
    u_exact = @(t) 3-3*t;
    x1_exact = @(t) (3./2)*t.^2-(1/2)*t.^3;
    x2_exact = @(t) 3*t - (3./2)*t.^2;
    
    % Plotting results
    newfig;
    plot(t_fwd, x1, 'r--', 'LineWidth', 1.5);
    hold on;
    plot(t_fwd, x1_exact(t_fwd), 'g', 'LineWidth', 1.5);
    hold on;
    plot(t_fwd, x2, 'b', 'LineWidth', 1.5);
    hold on;
    plot(t_fwd, x2_exact(t_fwd), 'm', 'LineWidth', 1.5);
    hold on;
    plot(t_lambda, u_new, 'k', 'LineWidth', 1.5);
    hold on; 
    plot(t_lambda, u_exact(t_lambda), 'm--', 'LineWidth', 1.5);
    hold on;
    xlabel('Time t');
    ylabel('Values');
    legend('x_1(t)','x1.exact(t)','x_2(t)','x2.exact(t)','u(t)','u.exact(t)');
    hold off;

my numerical and exact solutions are plotted below

[Torsten]

Maybe you could explain what you are trying to do mathematically in your code. I guess that most of the forum members have no experience with Pontryagin's maximum principle (like me).

Especially it would be interesting to know how you try to enforce the condition x1(1) = 1 in your code.

Answer 1

[shantanu]

Hi!

I understand you are attempting to implement the ‘PDEs FBSM method’ using ‘pontryagin principle’ through this MATLAB code and are then trying to compare the numerical solutions with the exact solutions.

I see there might be a few issues with your code:

The section where you are updating the temp variable values and are checking for convergence is commented out. Also go for a simpler convergence check condition if you feel the particular convergence technique you are using is never reached. Also in the interpolation parameter in state equation function, try to experiment with different techniques and choose the one that suits your usecase based on frequency of updates and compatibility with your algorithm.

It would help if you could share what is the expected output in this case so as to help figure where the present plot might be going wrong.

Here is one similar PDE I found: https://royalsocietypublishing.org/doi/10.1098/rsif.2021.0241 you may refer to this example.

Try to add breakpoints at places you are performing changes like calling state equation, costate equation, flips and see the values and accordingly try to debug would be the best way to resolve this error.

Hope this helps!