Newton's Method on MATLAB for Stationary Solutions for the Non-linear Klein-Gordon Equation

By interpreting the equation in this way, we can relate the dynamics described by the discrete Klein - Gordon equation to the behavior of DNA molecules within a biological system .

I am trying to represnt stationary points of the discrete Klein - Gordon equation. And the result are as follows.

Mathematica are different from MATLAB Any suggestion?

UPDATED CODE

% Parameters
numBases = 100; % Number of spatial points
omegaD = 0.2;  % Common parameter for the equation
% Preallocate the array for the function handles
equations = cell(numBases, 1);
% Initial guess for the solution
initialGuess = 0.01 * ones(numBases, 1);
% Parameter sets for kappa and beta
paramSets = [0.1, 0.05; 0.5, 0.05; 0.1, 0.2];
% Prepare figure for subplot
figure;
set(gcf, 'Position', [100, 100, 1200, 400]); % Set figure size

% Newton-Raphson method parameters
maxIterations = 1000;
tolerance = 1e-10;

for i = 1:size(paramSets, 1)
    kappa = paramSets(i, 1);
    beta = paramSets(i, 2);
    % Define the equations using a function
    for n = 2:numBases-1
        equations{n} = @(x) -kappa * (x(n+1) - 2 * x(n) + x(n-1)) - omegaD^2 * (x(n) - beta * x(n)^3);
    end
    % Boundary conditions with specified fixed values
    someFixedValue1 = 10; % Replace with actual value if needed
    someFixedValue2 = 10; % Replace with actual value if needed
    equations{1} = @(x) x(1) - someFixedValue1;
    equations{numBases} = @(x) x(numBases) - someFixedValue2;
    
    % Combine all equations into a single function
    F = @(x) cell2mat(cellfun(@(f) f(x), equations, 'UniformOutput', false));
    
    % Newton-Raphson iteration
    x_solution = initialGuess;
    for iter = 1:maxIterations
        % Evaluate function
        Fx = F(x_solution);
        
        % Evaluate Jacobian
        J = zeros(numBases, numBases);
        epsilon = 1e-6; % For numerical differentiation
        for j = 1:numBases
            x_temp = x_solution;
            x_temp(j) = x_temp(j) + epsilon;
            Fx_temp = F(x_temp);
            J(:, j) = (Fx_temp - Fx) / epsilon;
        end
        
        % Update solution
        delta = -J\Fx;
        x_solution = x_solution + delta;
        
        % Check convergence
        if norm(delta, inf) < tolerance
            fprintf('Convergence achieved after %d iterations.\n', iter);
            break;
        end
    end
    
    if iter == maxIterations
        error('Maximum iterations reached without convergence.');
    end
    
    % Plot the solution in a subplot
    subplot(1, 3, i);
    plot(x_solution, 'o-', 'LineWidth', 2);
    grid on;
    xlabel('n', 'FontSize', 12);
    ylabel('x[n]', 'FontSize', 12);
    title(sprintf('\\kappa = %.2f, \\beta = %.2f', kappa, beta), 'FontSize', 14);
end

Convergence achieved after 392 iterations.
Convergence achieved after 424 iterations.
Convergence achieved after 63 iterations.

% Improve overall aesthetics
sgtitle('Stationary States for Different \kappa and \beta Values', 'FontSize', 16); % Super title for the figure

Mathematica plots

https://community.wolfram.com/groups/-/m/t/3137876

[VBBV]

try changing this parameter and check

omegaD = 0.2;  % Common parameter for the equation

What is common parameter for the equation represent ?

[Athanasios Paraskevopoulos]

representing the frequency term, speculative and low to reflect slow dynamics

[Athanasios Paraskevopoulos]

Maybe they have different precision

Answer 1

[Torsten]

All you see in the graphs are floating point errors.

Your system of equations has solution x(i) = 0 for 1 <= i <= numBases and all three parameter constellations.

[Athanasios Paraskevopoulos]

@Torsten the code I tried to translate to Matlab code is the following https://community.wolfram.com/groups/-/m/t/3137876. What should I change to avoid floating point errors?

[Torsten]

Isn't your system of equations the same in MATLAB and in MATHEMATICA ? My guess is that you see floating point trash in both software environments.

[Athanasios Paraskevopoulos]

@Torsten Yes, I solved the same system! Then I have to wonder why they chose my code for their editorial column if it's trash

[Torsten]

Maybe they didn't recognize that x = 0 solves your system of equations.

At least it will be very difficult to prove with numerical methods that there is another solution in the order of 1e-19 besides 0.

Can't you argue from a physical perspective that a boundary condition of 0 at both sides will result in x = 0 ?

[Athanasios Paraskevopoulos]

@Torsten No I can not! if I change the boundary condition could I receive more realistic results?

[Torsten]

Without the nonlinear term omegaD^2 * (x(n) - beta * x(n)^3), x = 0 is the solution of your equations - thus very realistic. It's also the stationary solution of the heat equation where both ends are held at temperature 0 - which will lead to a temperature 0 in the complete region - physical reasoning.

I don't want to exclude the possibility that solutions apart from 0 could exist, but - as said - it will be hard to prove if they are in the order of magnitude as in your example cases.

Sorry, but I have no experience with the Klein-Gordon equation. Of course, solutions will change if you set one or both of the boundary conditions to a number different from 0. But I couldn't interprete if the results thus obtained are realistic or not and if they have a physical interpretation.

[Athanasios Paraskevopoulos]

Thank you very much for your time

[Torsten]

If you start with

initialGuess = 5 * ones(numBases, 1);

you definitly get solutions different from 0.

[Athanasios Paraskevopoulos]

The results are more realistic now. Like a snapshot of DNA https://www.mathworks.com/matlabcentral/discussions/general/848831-exploring-dna-oscillations-with-the-discrete-klein-gordon-framework

[Athanasios Paraskevopoulos]

@Torsten I tried Newton's Method on both softwares but the results are so different

[Torsten]

Did you test whether the solutions obtained satisfy your equations ?

Because your convergence check is only on norm(x-xold) < tol, not on norm(f(x)) < tol.

I did this for you, and the solutions satisfy your equations. So from MATLAB side, you got a solution (that might depend on the initial values you prescribe).

[Athanasios Paraskevopoulos]

Is this the code that you mean?

% Parameters
numBases = 100; % Number of spatial points
omegaD = 0.2;  % Common parameter for the equation
% Preallocate the array for the function handles
equations = cell(numBases, 1);
% Initial guess for the solution
initialGuess = 0.01 * ones(numBases, 1);
% Parameter sets for kappa and beta
paramSets = [0.1, 0.05; 0.5, 0.05; 0.1, 0.2];
% Prepare figure for subplot
figure;
set(gcf, 'Position', [100, 100, 1200, 400]); % Set figure size

% Newton-Raphson method parameters
maxIterations = 1000;
tolerance = 1e-10;

for i = 1:size(paramSets, 1)
    kappa = paramSets(i, 1);
    beta = paramSets(i, 2);
    % Define the equations using a function
    for n = 2:numBases-1
        equations{n} = @(x) -kappa * (x(n+1) - 2 * x(n) + x(n-1)) - omegaD^2 * (x(n) - beta * x(n)^3);
    end
    % Boundary conditions with specified fixed values
    someFixedValue1 = 10; % Replace with actual value if needed
    someFixedValue2 = 10; % Replace with actual value if needed
    equations{1} = @(x) x(1) - someFixedValue1;
    equations{numBases} = @(x) x(numBases) - someFixedValue2;
    
    % Combine all equations into a single function
    F = @(x) cell2mat(cellfun(@(f) f(x), equations, 'UniformOutput', false));
    
    % Newton-Raphson iteration
    x_solution = initialGuess;
    for iter = 1:maxIterations
        % Evaluate function
        Fx = F(x_solution);
        
        % Evaluate Jacobian
        J = zeros(numBases, numBases);
        epsilon = 1e-6; % For numerical differentiation
        for j = 1:numBases
            x_temp = x_solution;
            x_temp(j) = x_temp(j) + epsilon;
            Fx_temp = F(x_temp);
            J(:, j) = (Fx_temp - Fx) / epsilon;
        end
        
        % Update solution
        delta = -J\Fx;
        x_solution = x_solution + delta;
        
        % Check convergence based on solution change and function value
        if norm(delta, inf) < tolerance && norm(Fx, inf) < tolerance
            fprintf('Convergence achieved after %d iterations.\n', iter);
            break;
        end
    end
    
    if iter == maxIterations
        error('Maximum iterations reached without convergence.');
    end
    
    % Plot the solution in a subplot
    subplot(1, 3, i);
    plot(x_solution, 'o-', 'LineWidth', 2);
    grid on;
    xlabel('n', 'FontSize', 12);
    ylabel('x[n]', 'FontSize', 12);
    title(sprintf('\\kappa = %.2f, \\beta = %.2f', kappa, beta), 'FontSize', 14);
end

Convergence achieved after 392 iterations.
Convergence achieved after 424 iterations.
Convergence achieved after 63 iterations.

% Improve overall aesthetics
sgtitle('Stationary States for Different \kappa and \beta Values', 'FontSize', 16); % Super title for the figure

[Torsten]

if norm(delta, inf) < tolerance && norm(Fx, inf) < tolerance

Yes, that's what I meant.

And "fsolve" gives different results ?

[Athanasios Paraskevopoulos]

@Torsten Why sometimes MatLab and Mathematica lead to different results is intriguing me. If both are based on numerical methods very well tested, both should arrive to the same results. As you can see here the plots are totally different https://community.wolfram.com/groups/-/m/t/3137876

[Torsten]

If you use the code below, you will see that you get similar results as in MATHEMATICA, but that the residual is too high to be converged. You should check norm(F(solution)) or a possible error flag from the solver in MATHEMATICA to see whether you can trust the results.

% Parameters
numBases = 100; % Number of spatial points
omegaD = 0.2;  % Common parameter for the equation
% Preallocate the array for the function handles
equations = cell(numBases, 1);
% Initial guess for the solution
initialGuess = 0.01 * ones(numBases, 1);
% Parameter sets for kappa and beta
paramSets = [0.1, 0.05; 0.5, 0.05; 0.1, 0.2];
% Prepare figure for subplot
figure;
set(gcf, 'Position', [100, 100, 1200, 400]); % Set figure size

% Newton-Raphson method parameters
maxIterations = 1000;
tolerance = 1e-10;

for i = 1:size(paramSets, 1)
    kappa = paramSets(i, 1);
    beta = paramSets(i, 2);
    % Define the equations using a function
    for n = 2:numBases-1
        equations{n} = @(x) -kappa * (x(n+1) - 2 * x(n) + x(n-1)) - omegaD^2 * (x(n) - beta * x(n)^3);
    end
    % Boundary conditions with specified fixed values
    someFixedValue1 = 10; % Replace with actual value if needed
    someFixedValue2 = 10; % Replace with actual value if needed
    equations{1} = @(x) x(1) - someFixedValue1;
    equations{numBases} = @(x) x(numBases) - someFixedValue2;
    
    % Combine all equations into a single function
    F = @(x) cell2mat(cellfun(@(f) f(x), equations, 'UniformOutput', false));
    
    x_solution = fsolve(F,initialGuess);
    norm(F(x_solution))
    
    % Plot the solution in a subplot
    subplot(1, 3, i);
    plot(x_solution, 'o-', 'LineWidth', 2);
    grid on;
    xlabel('n', 'FontSize', 12);
    ylabel('x[n]', 'FontSize', 12);
    title(sprintf('\\kappa = %.2f, \\beta = %.2f', kappa, beta), 'FontSize', 14);
end

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 1.000000e+04.

ans = 0.0966

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 1.000000e+04.

ans = 0.1420

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 1.000000e+04.

ans = 0.1437

% Improve overall aesthetics
sgtitle('Stationary States for Different \kappa and \beta Values', 'FontSize', 16); % Super title for the figure

[Athanasios Paraskevopoulos]

Because I am begginer to Matlab, when I use fsolve

% Solve the system of equations using fsolve
x_solution = fsolve(F, initialGuess);
norm(F(x_solution));

I replace all the following code and is more accurate right?

% Newton-Raphson iteration
for iter = 1:maxIterations
    % Evaluate function
    Fx = F(x_solution);
    
    % Evaluate Jacobian
    J = zeros(numBases, numBases);
    epsilon = 1e-6; % For numerical differentiation
    for j = 1:numBases
        x_temp = x_solution;
        x_temp(j) = x_temp(j) + epsilon;
        Fx_temp = F(x_temp);
        J(:, j) = (Fx_temp - Fx) / epsilon;
    end
    
    % Update solution
    delta = -J\Fx;
    x_solution = x_solution + delta;
    
    % Check convergence
    if norm(delta, inf) < tolerance && norm(Fx, inf) < tolerance
        fprintf('Convergence achieved after %d iterations.\n', iter);
        break;
    end
end

if iter == maxIterations
    error('Maximum iterations reached without convergence.');
end

[Torsten]

Because I am begginer to Matlab, when I use fsolve I replace all the following code

Yes.

and is more accurate right?

Obviously not in your case since it doesn't converge for the initialGuess you supply. Only if you set

initialGuess = 10 * ones(numBases, 1);

e.g., it converges to the solution you get with your code of Newton's method (except for the first case where at least three solutions seem to exist: 0, your solution and the solution you get with "fsolve").

So it's hard to interprete the results if you cannot apply physical reasoning to sort out solutions that make no sense.

Answer 2

[Athanasios Paraskevopoulos]

@Torsten I read the fsolve Algorithms and here's the modified MATLAB script using the levenberg-marquardt algorithm for the fsolve function. Now I have the results that I wanted:

% Parameters
numBases = 100; % Number of spatial points
omegaD = 0.2;  % Common parameter for the equation
% Preallocate the array for the function handles
equations = cell(numBases, 1);
% Initial guess for the solution
initialGuess = 0.01 * ones(numBases, 1);
% Parameter sets for kappa and beta
paramSets = [0.1, 0.05; 0.5, 0.05; 0.1, 0.2];
% Prepare figure for subplot
figure;
set(gcf, 'Position', [100, 100, 1200, 400]); % Set figure size

% Newton-Raphson method parameters
maxIterations = 1000;
tolerance = 1e-10;

% Set options for fsolve to use the 'levenberg-marquardt' algorithm
options = optimoptions('fsolve', 'Algorithm', 'levenberg-marquardt', 'MaxIterations', maxIterations, 'FunctionTolerance', tolerance);

for i = 1:size(paramSets, 1)
    kappa = paramSets(i, 1);
    beta = paramSets(i, 2);
    % Define the equations using a function
    for n = 2:numBases-1
        equations{n} = @(x) -kappa * (x(n+1) - 2 * x(n) + x(n-1)) - omegaD^2 * (x(n) - beta * x(n)^3);
    end
    % Boundary conditions with specified fixed values
    someFixedValue1 = 10; % Replace with actual value if needed
    someFixedValue2 = 10; % Replace with actual value if needed
    equations{1} = @(x) x(1) - someFixedValue1;
    equations{numBases} = @(x) x(numBases) - someFixedValue2;
    
    % Combine all equations into a single function
    F = @(x) cell2mat(cellfun(@(f) f(x), equations, 'UniformOutput', false));
    
    % Solve the system of equations using fsolve with the specified options
    x_solution = fsolve(F, initialGuess, options);
    norm(F(x_solution))
    
    % Plot the solution in a subplot
    subplot(1, 3, i);
    plot(x_solution, 'o-', 'LineWidth', 2);
    grid on;
    xlabel('n', 'FontSize', 12);
    ylabel('x[n]', 'FontSize', 12);
    title(sprintf('\\kappa = %.2f, \\beta = %.2f', kappa, beta), 'FontSize', 14);
end

No solution found.

fsolve stopped because the last step was ineffective. However, the vector of function
values is not near zero, as measured by the value of the function tolerance. 

ans = 0.0569

No solution found.

fsolve stopped because the last step was ineffective. However, the vector of function
values is not near zero, as measured by the value of the function tolerance. 

ans = 0.1339

No solution found.

fsolve stopped because the relative size of the current step is less than the
value of the step size tolerance, but the vector of function values
is not near zero as measured by the value of the function tolerance.

ans = 0.0306

% Improve overall aesthetics
sgtitle('Stationary States for Different \kappa and \beta Values', 'FontSize', 16); % Super title for the figure