Need Help with the code. Why I am getting the zero value for V., G_V and K_V?

The system of PDE which i have is given by

where

Also the boundary conditions are

I have a written a code

% Define parameters
A_star = 60;   % Upper limit for a
T = 1.0;        % Time integration limit
L = 10.0;       % Spatial domain limit
d = 1.0;        % Diffusion coefficient
lambda_tau = 1.0;
U0 = 1.0;

% Discretization
Ny = 100; Na = 100;
dy = L / Ny; da = A_star / Na;
y = linspace(-L, L, Ny);  % Spatial grid
a = linspace(0, A_star, Na);  % Age grid

% Initialize U and V
U = U0 * ones(Ny, 1);  % Initial value for U
V = zeros(Na, Ny);      % Initial value for V

% Define kernel functions
K = @(x) exp(-x.^2);  % Example kernel function
G = @(tau, s) exp(-tau.^2 - s.^2);  % Example G function
eta = @(a) exp(-a);  % Example eta function

% Convolution function
function result = convolve_K(f, K, dy)
    N = length(f);
    result = zeros(N, 1);
    for i = 1:N
        for j = 1:N
            result(i) = result(i) + K(abs(i - j) * dy) * f(j) * dy;
        end
    end
end

% G convolution for V
function result = G_convolution(V, G, y, a, lambda_tau, dy)
    [Na, Ny] = size(V);
    result = zeros(Na, Ny);
    for i = 1:Na
        for j = 1:Ny
            sum_val = 0;
            for tau = 1:Ny
                for s = 1:Ny
                    y_idx = j - round(lambda_tau * tau * dy);
                    if y_idx > 0 && y_idx <= Ny
                        sum_val = sum_val + G(tau * dy, s * dy) * V(i, y_idx) * dy^2;
                    end
                end
            end
            result(i, j) = sum_val;
        end
    end
end

% Time-stepping parameters
dt = 0.01;  % Time step
n_steps = 100;  % Number of time steps

% Time stepping loop
for t = 1:n_steps
    % Convolution of U with K
    K_U = convolve_K(U, K, dy);

    % Compute integral term for U update
    G_V = G_convolution(V, G, y, a, lambda_tau, dy);
    integral_term = zeros(Ny, 1);
    for j = 1:Ny
        integral_term(j) = trapz(a, eta(a(:)) .* G_V(:, j));
    end

    % Update U
    S = ones(Ny, 1);  % Example S(y)
    dU_dy = d * (K_U - U) - S .* integral_term;
    U = U + dt * dU_dy;

    % Convolution of V with K
    K_V = zeros(size(V));
    for i = 1:Na
        K_V(i, :) = convolve_K(V(i, :)', K, dy)';
    end

    % Update V
    dV_dy = d * (K_V - V);
    dV_da = diff(V, 1, 1) / da;
    dV_da = [zeros(1, Ny); dV_da];  % Boundary condition for a = 0
    V = V + dt * (dV_dy - dV_da);

    % Apply boundary condition for V(0, y)
    V(1, :) = U' .* trapz(a, eta(a) .* G_V, 1);
end

% Plot results
figure;
subplot(2, 1, 1);
plot(y, U);
title('U(y)');
xlabel('y'); ylabel('U');

subplot(2, 1, 2);
surf(y, a, V,'EdgeColor','none');
axis([-10 10 0 60 0 1])
xlabel('space');
ylabel('Age');
zlabel('Educated Susceptibles');
colormap jet    % change color map
colorbar   

This code works perfectly for $U$, but the value of V, G_V, K_V is coming out to be zero everywhere. What mistake did I made? Is there any simpler way to write code for my model. Please help me with this.

Answer 1

[Divyajyoti Nayak]

Hi @Nikhil,

The reason for why the values of matrix “K_V” is coming out to be zero is because the values are being initialized as zero, then multiplied with other variables which gives zero and added to initial “result”. Hence, the new value of “result” is the same as the old value of “result” which was also zero. Here’s the part of the code that I am referring to:

% Convolution function
function result = convolve_K(f, K, dy)
N = length(f);
result = zeros(N, 1);
    for i = 1:N
        for j = 1:N
            result(i) = result(i) + K(abs(i - j) * dy) * ...
                f(j) ... f(j) = K_V(i,j) = 0 (Initialized value)
                * dy; %Hence result(i) = 0 (always)
        end
    end
end

Similarly, the values of “V” and “G_V” are also zero.