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

Question

Athanasios Paraskevopoulos 2024-3-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2095571-newton-s-method-on-matlab-for-stationary-solutions-for-the-non-linear-klein-gordon-equation

编辑： Torsten 2024-3-20

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

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

Athanasios Paraskevopoulos 2024-3-18

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

Athanasios Paraskevopoulos 2024-3-18

Maybe they have different precision

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-3-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2095571-newton-s-method-on-matlab-for-stationary-solutions-for-the-non-linear-klein-gordon-equation#answer_1427226

编辑：Torsten 2024-3-18

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.

17 个评论
显示 15更早的评论隐藏 15更早的评论

Athanasios Paraskevopoulos 2024-3-20

在 MATLAB Online 中打开

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

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 2024-3-20

在 MATLAB Online 中打开

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 2024-3-20

在 MATLAB Online 中打开

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 2024-3-20

编辑：Torsten 2024-3-20

在 MATLAB Online 中打开

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 2024-3-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2095571-newton-s-method-on-matlab-for-stationary-solutions-for-the-non-linear-klein-gordon-equation#answer_1428531

在 MATLAB Online 中打开

@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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

采纳的回答

17 个评论
显示 15更早的评论隐藏 15更早的评论

更多回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

3 个评论 显示 1更早的评论隐藏 1更早的评论

采纳的回答

17 个评论 显示 15更早的评论隐藏 15更早的评论

更多回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

3 个评论
显示 1更早的评论隐藏 1更早的评论

17 个评论
显示 15更早的评论隐藏 15更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论