Why numerical solutions is not consistent with the exact solutions?

Why numerical solutions is not consistent with the exact solutions? Also the error is too large. How to modify the code to get consistent results? See the attached image for the equation of motions and their exact solutions.

clear; clc;
tic;  
c1 = 0; % real constant
c2 = 0; % real constant
a = 1; % real constant
b=0.1; % two values b=0.1, b=0.01
lambda = -1./(2*(1 + 1i)); % complex or imaginary
N = 1000;
X = linspace(-10, 10, N);
T = linspace(0, 10, N);

% Calculate initial conditions
[x, ~] = meshgrid(X, 0); % t = 0 for initial conditions
%psi_init = A * exp(-1i * lambda * x);
%phi_init = B * exp(1i * lambda * x);
nq_init = -a*(b.^5*(1+(0 + c2).^2).^2 - 8*a*(0+c2).*(x+c1) - 4*a*b.^4*(0+c2).*(1+(0+c2).^2).*(x+c1) - 4*a*b.^2*(0+c2).*(x+c1).*(3+(0+c2).^2 + a.^2*(x+c1).^2) + b*(-3+(0+c2).^4 + 6*a.^2*(x+c1).^2 + a.^4*(x+c1).^4 + 2*(0+c2).^2.*(3 + a.^2*(x+c1).^2)) + 2*b.^3*(-1 + (0+c2).^4 + a.^2*(x+c1).^2 + (0+c2).^2*(4 + 3*a.^2*(x+c1).^2)));
dq_init = 1 + (0 + c2).^2 + b.^2 + (b*(0 + c2) + a*(x + c1)).^2;
q_init = (nq_init ./ (dq_init .^ 2));
ns_init = -(-3 + (0 + c2).^2 + b.^2*(-3 + (0 + c2).^2) + 4*1i*a*(x+c1) - 2*a*b*(0 + c2).*(x+c1) + a.^2*(x+c1).^2).*(exp(1i*(a*x + b*0)));
s_init = ns_init ./ dq_init;

initial_conditions = [q_init(:); s_init(:)];

t_span = T; 

% Solve the differential equation numerically
options = odeset('RelTol', 1e-4, 'AbsTol', 1e-8);
[t, u] = ode45(@(t, u) dydt(t, u, N), t_span, initial_conditions, options);

% Extract q and s from the solution
q = reshape(u(:, 1:N), [length(t), N]);
s = reshape(u(:, N+1:end), [length(t), N]);

% Calculate the exact solution over the mesh grid
[x, t] = meshgrid(X, T);
nq = -a*(b.^5*(1+(t + c2).^2).^2 - 8*a*(t+c2).*(x+c1) - 4*a*b.^4*(t+c2).*(1+(t+c2).^2).*(x+c1) - 4*a*b.^2*(t+c2).*(x+c1).*(3+(t+c2).^2 + a.^2*(x+c1).^2) + b*(-3+(t+c2).^4 + 6*a.^2*(x+c1).^2 + a.^4*(x+c1).^4 + 2*(t+c2).^2.*(3 + a.^2*(x+c1).^2)) + 2*b.^3*(-1 + (t+c2).^4 + a.^2*(x+c1).^2 + (t+c2).^2*(4 + 3*a.^2*(x+c1).^2)));
dq = 1 + (t + c2).^2 + b.^2 + (b*(t + c2) + a*(x + c1)).^2;
ns = -(-3 + (t + c2).^2 + b.^2*(-3 + (t + c2).^2) + 4*1i*a*(x+c1) - 2*a*b*(t + c2).*(x+c1) + a.^2*(x+c1).^2).*(exp(1i*(a*x + b*0)));
exact_q = (nq ./ (dq .^ 2));
exact_s = ns ./ dq;

% Compute the error
error_q = abs(q' - exact_q);
error_s = abs(s' - exact_s);

% Plotting the numerical solution
figure;
subplot(2, 1, 1);
surf(t, x, abs(q'));
title('Numerical Solution for |q_{n+1} - q_n|');
xlabel('t');
ylabel('x');
zlabel('|q_{n+1} - q_n|');
shading interp;

subplot(2, 1, 2);
surf(t, x, abs(s'));
title('Numerical Solution for |s_{n+1} - s_n|');
xlabel('t');
ylabel('x');
zlabel('|s_{n+1} - s_n|');
shading interp;

% 3D surface plot for exact solution
figure;
subplot(2, 1, 1);
surf(t, x, abs(exact_q'));
title('Exact Solution for |q_{n+1} - q_n|');
xlabel('t');
ylabel('x');
zlabel('|q_{n+1} - q_n|');
shading interp;

subplot(2, 1, 2);
surf(t, x, abs(exact_s'));
title('Exact Solution for |s_{n+1} - s_n|');
xlabel('t');
ylabel('x');
zlabel('|s_{n+1} - s_n|');
shading interp;

% Plot the error
figure;
subplot(2, 1, 1);
surf(t, x, error_q');
title('Error for |q_{n+1} - q_n|');
xlabel('t');
ylabel('x');
zlabel('Error');
shading interp;

subplot(2, 1, 2);
surf(t, x, error_s');
title('Error for |s_{n+1} - s_n|');
xlabel('t');
ylabel('x');
zlabel('Error');
shading interp;

% 2D error plot at specific times
figure;
times_to_plot = [0, 5, 10]; % Example times to plot error
for i = 1:length(times_to_plot)
    time_idx = find(t_span >= times_to_plot(i), 1);
    subplot(length(times_to_plot), 1, i);
    plot(X, error_q(:, time_idx), 'r', X, error_s(:, time_idx), 'b');
    title(['Error at t = ', num2str(times_to_plot(i))]);
    xlabel('x');
    ylabel('Error');
    legend('|q_{n+1} - q_n| Error', '|s_{n+1} - s_n| Error');
end

% Generate error table for fixed times from 0 to 10 in increments of 0.5
fixed_times = 0:0.5:10;
error_table_q = zeros(length(fixed_times), 1);
error_table_s = zeros(length(fixed_times), 1);

for i = 1:length(fixed_times)
    time_idx = find(t_span >= fixed_times(i), 1);
    error_table_q(i) = mean(error_q(:, time_idx));
    error_table_s(i) = mean(error_s(:, time_idx));
end

% Display the error table
error_table = table(fixed_times', error_table_q, error_table_s, ...
                    'VariableNames', {'Time', 'Error_q', 'Error_s'});
disp(error_table);

    Time    Error_q    Error_s
    ____    _______    _______

      0      1.0826    1.4923 
    0.5      1.1705    1.3553 
      1      1.2549    1.5695 
    1.5      1.3289    1.8461 
      2      1.3825    1.9677 
    2.5      1.4034    1.8828 
      3      1.3782     1.625 
    3.5      1.2962    1.4018 
      4      1.1567    1.4904 
    4.5     0.95671    1.6815 
      5     0.54434    1.7139 
    5.5      0.4193    1.6455 
      6      0.5698    1.4517 
    6.5     0.70902    1.3707 
      7      0.7934    1.5314 
    7.5     0.83008    1.7021 
      8     0.83204    1.7524 
    8.5     0.81079    1.6522 
      9     0.77526    1.4462 
    9.5     0.73188    1.3023 
     10     0.68601    1.4173 

toc;

Elapsed time is 1.211971 seconds.

function dydt = dydt(~, u, N)
    dydt = zeros(2*N, 1);
    q = u(1:N);
    s = u(N+1:2*N);

    % Handle boundary conditions
    for n = 3:N-2
        dydt(n) = 1/2 * s(n) * conj(s(n)) - 1/2 * s(n+1) * conj(s(n+1));
        dydt(N+n) = 1/2 * (q(n+1) - q(n)) * (s(n+1) - s(n));
    end
end

[Torsten]

What are the underlying partial differential equations for q and s you try to solve ?

What are the boundary conditions you try to impose ?

The loop

for n = 3:N-2
        dydt(n) = 1/2 * s(n) * conj(s(n)) - 1/2 * s(n+1) * conj(s(n+1));
        dydt(N+n) = 1/2 * (q(n+1) - q(n)) * (s(n+1) - s(n));
end

does not assign time derivatives to q(1),q(2),q(N-1),q(N),s(1),s(2),s(N-1) and s(N). Is this intentionally or just an error ?

[A]

PDEs are

q is real function, and s is a complex function, represents complex conjugate of s.

The system has natural boundary conditions:

i) I used initial condition rather boundary conditions.

ii) Means "dydt(n)" should be "dydtq(n)", and "dydt(N+n)" should be "dydts(N+n)"?

[Torsten]

I don't understand your comments i) and ii).

Why you plot the difference of the solution in neighbouring cells ?

Since the solution seems to be symmetric, couldn't you solve it from 0 to x with boundary condition s_x = q_x = 0 at x = 0 ?

Why is there no dt or dx in your equations ? Is it natural to choose both as 1 ?

[A]

@Torsten

1. Why do you plot the difference of the solution in neighbouring cells?

We plot the differences because the governing equations are naturally written in terms of these differences:

This structure comes from the nature of the system (see the attached image at the top of the page) being studied.

Solving for these differences rather than the full field values $q_n(t)$ and $s_n(t)$ allows us to directly follow the dynamics as they are defined in the model. These terms capture the local interactions between adjacent points on the spatial lattice, much like fluxes or gradients in finite difference methods for partial differential equations.

2. Since the solution seems symmetric, couldn't you solve it from with boundary conditions ?

In theory, yes. If the solution is known to be symmetric about , we could reduce the domain to and impose Neumann boundary conditions:

However, in practice, we solve the problem on the symmetric domain to preserve full generality and avoid introducing artificial constraints or boundary effects. This is particularly important when:

i. The initial condition is not perfectly symmetric.

ii. We later introduce perturbations or external influences.

iii. We want to verify symmetry emerges naturally from the dynamics.

3. Why is there no or in your equations? Is it natural to choose both as 1?

This model is semi-discrete: it is discrete in space (indexed by n) and continuous in time (t). Therefore:

i. There is no explicit spatial step in the equations. It is implicitly assumed that the lattice spacing is uniform and normalized to 1. For example:

ii. Time is treated continuously, and we use an ODE solver (such as ode15s) to integrate over time. Hence, no explicit appears in the equations.

Thus, yes — it is natural to assume and work in dimensionless units. If needed, physical dimensions can be restored through appropriate rescaling.

[Torsten]

In function "dydt", you supply the time derivatives for q_x and s_x.

Consequently, your solution variables are q_x and s_x, not q and s.

Thus

q = u(1:N);
s = u(N+1:2*N);  

is wrong.

This should somehow be

q_x = u(1:N-1);
s_x = u(N:2*N-2);   

And initial conditions also have to be supplied for q_x and s_x, not q and s.

And be careful with the ' - operator. For complex vectors, ' means conjugate-transpose, not only transpose:

z = [3 + 1i*5, -2 - 1i*6];
z.'

ans = 
   3.0000 + 5.0000i
  -2.0000 - 6.0000i

z'

ans = 
   3.0000 - 5.0000i
  -2.0000 + 6.0000i

[A]

@Torsten I made the corrections as you suggested. But it's not working properly and I couldn't fix it yet. See the attached files.