Non-homogenous spatial and temporal grids

Good day family. I am solving a system of 2 PDEs using the pdepe package (The algorithm is in the matlab documentation: https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html ). The parameters of the system depend on some observations (data) and the spatial axis as well as shown in my code below:

% beggining of the code
function non_constant_spatial_and_temporal_axis
% pdebc function -- Neumann BCs
    function [pl, ql, pr, qr] = pdebc(~, ~, ~, ~, ~)
        pl = [0; 0];
        ql = [1; 1];
        pr = [0; 0];
        qr = [1;1];
    end
% pdeic function
    function u0 = pdeic(~)
        u0 = [10;5];
    end
% pdefun function
    function [c, f, s] = pdefun(x, ~, u, dudx)
        D1 = 0.024;
        D2 = 0.0170;
        c = [1; 1];
        f = [D1; D2].*dudx;
        
        t1=[1,0,3,4,5,3,2,3,2,3,3,2,3,4,5,6,7,5,4,3,3,2,2,1,3,2,4,5]; % This data is the observations which represent the temporal grid
        x_interpolate = linspace(0, 10, length(t1));     % Interpolation of the data so that it can match with the spatial grid
        t1_interpolated = interp1(x_interpolate,t1,x,'linear');
        alpha1 = 0.5*(1-sin(2*t1_interpolated));  % The parameters alpha1 and alpha2 depends on the observations in t1
        % I want to define alpha1 in such a way that it is a function of bth
        % time and space x such that alpha1(t,x)=
        % 0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x)) so that I can visualize
        % the oscillations on both the time and the space axis at the same
        % time. 
        alpha2 = 0.5*(1-sin(2*t1_interpolated));   % The parameter alpha2 depends on the observations in t1
        % alpha2(t,x)=0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x));
        s = [alpha1.*u(2).*(u(2).^2 -1); -alpha2.* u(1).*(u(2)-1)];
    end
% Main script for running the code
clear all; close all; clc;
m = 0;

t = linspace(0, 5, 50); % time variable
x = linspace(0, 10,50); % space varable
sol = pdepe(m, @pdefun, @pdeic, @pdebc, x, t);
% Extract the solutions
u1 = sol(:, :, 1);
u2 = sol(:, :, 2);
% Plots
figure(1)
surf(x, t, u1)
title('u_1(x,t)')
xlabel('Time t')
ylabel('Postion x')
figure(2)
surf(x, t, u2)
title('u_2(x,t)')
xlabel('Time t')
ylabel('Position x')
end
%end of code

The code is running well but its not giving me oscillations on the spatial axis, like what I am getting on the temporal axis. I wish to be able to observe the spatial variation using a sinusoidal function like the sin function and I also view the oscillations on the spatial and temporal axis. Note that the oscillations on the temporal axis are due to the data while the ones on the spatial axis are due to the use of the sinusoidal function. May you kindly assist me resolve this. I have tried substituting alpha1 with alpha1=0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x)); but then I observe that the term (1-sin(2*x)) reduces the value of alpha1 and the effect is only felt on the temporal axis and not in the spatial axis. What can I do to have the effect felt on the spatial axis as well?

Answer 1

[Sanju]

Hi @Mr. P,

In the “pdefun” function, you want to define “alpha1” and “alpha2” as functions of both time and space. You can achieve this by modifying the definition of “alpha1” and “alpha2” as follows:

alpha1 = 0.5*(1-sin(2*t1_interpolated)).*(1-sin(2*x));
alpha2 = 0.5*(1-sin(2*t1_interpolated)).*(1-sin(2*x));

This will make “alpha1” and “alpha2” functions of both t and x, and you can visualize the oscillations on both the time and space axes.

To observe oscillations on the spatial axis, you can modify the definition of “alpha1” and “alpha2” as follows:

alpha1 = 0.5*(1-sin(2*t1_interpolated)).*(1-sin(2*pi*x/10));
alpha2 = 0.5*(1-sin(2*t1_interpolated)).*(1-sin(2*pi*x/10));

This will make “alpha1” and “alpha2” functions of both t and x, and you can visualize the oscillations on both the time and space axes. The 2*pi*x/10 term in the definition of “alpha1” and “alpha2” introduces a sinusoidal variation in the spatial direction.

The x/10 term is used to normalize the spatial axis to the range [0,1].

The difference between using x and x/10 is that x represents the spatial axis in the original units, while x/10 represents the spatial axis normalized to the range [0,1]. Therefore, x/10 is used to normalize the spatial axis to the required range.

Here’s is the updated code,

% beggining of the code
function non_constant_spatial_and_temporal_axis
    % pdebc function -- Neumann BCs
function [pl, ql, pr, qr] = pdebc(~, ~, ~, ~, ~)
    pl = [0; 0];
    ql = [1; 1];
    pr = [0; 0];
    qr = [1;1];
end
% pdeic function
function u0 = pdeic(~)
    u0 = [10;5];
end
 % pdefun function
function [c, f, s] = pdefun(x, ~, u, dudx)
    D1 = 0.024;
    D2 = 0.0170;
    c = [1; 1];
    f = [D1; D2].*dudx;
    
    t1=[1,0,3,4,5,3,2,3,2,3,3,2,3,4,5,6,7,5,4,3,3,2,2,1,3,2,4,5]; % This data is the observations which represent the temporal grid
    x_interpolate = linspace(0, 10, length(t1));     % Interpolation of the data so that it can match with the spatial grid
    t1_interpolated = interp1(x_interpolate,t1,x,'linear');


   alpha1 = 0.5*(1-sin(2*t1_interpolated)).*(1-sin(2*pi*x/10));
   alpha2 = 0.5*(1-sin(2*t1_interpolated)).*(1-sin(2*pi*x/10));

    %alpha1 = @(x,t,u,dudx) 0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x));
   % alpha1 = 0.5*(1-sin(2*t1_interpolated));  % The parameters alpha1 and alpha2 depends on the observations in t1
    % I want to define alpha1 in such a way that it is a function of bth
    % time and space x such that alpha1(t,x)=
    % 0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x)) so that I can visualize
    % the oscillations on both the time and the space axis at the same
    % time. 
   % alpha2 = 0.5*(1-sin(2*t1_interpolated));   % The parameter alpha2 depends on the observations in t1
    % alpha2(t,x)=0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x));
    s = [alpha1.*u(2).*(u(2).^2 -1); -alpha2.* u(1).*(u(2)-1)];
end
    % Main script for running the code
    clear all; close all; clc;
    m = 0;
     
    t = linspace(0, 5, 50); % time variable
    x = linspace(0, 10,50); % space varable
    sol = pdepe(m, @pdefun, @pdeic, @pdebc, x, t);
    % Extract the solutions
    u1 = sol(:, :, 1);
    u2 = sol(:, :, 2);
    % Plots
    figure(1)
    surf(x, t, u1)
    title('u_1(x,t)')
    xlabel('Time t')
    ylabel('Postion x')
    figure(2)
    surf(x, t, u2)
    title('u_2(x,t)')
    xlabel('Time t')
    ylabel('Position x')
end
%end of code

Hope this Helps!

Thanks.

[Mr. P]

Dear @Sanju Thank you very much for your response. However, my problem is to view the oscillations on both time and space at the same time, without having to view the other axis as constant. Any ideas?

[Sanju]

One way to view oscillations on both time and space simultaneously is by creating a 3D plot. You can use the meshgrid function to generate a grid of points in both time and space, and then plot the oscillations as the height of the surface.

Here's an example code snippet that demonstrates this approach,

% Define the time and space vectors
t = linspace(0, 10, 100); % Time vector
x = linspace(0, 1, 100); % Space vector

% Create a grid of points in time and space
[tt, xx] = meshgrid(t, x);

% Compute the oscillations at each point in the grid
oscillations = sin(2*pi*tt) .* sin(2*pi*xx);

% Create a 3D plot
figure;
surf(tt, xx, oscillations);
xlabel('Time');
ylabel('Space');
zlabel('Oscillations');

This code will generate a 3D plot where the x-axis represents time, the y-axis represents space, and the z-axis represents the oscillations. You can rotate the plot to view the oscillations from different angles.