Error for nonconstant Neumann boundary condition when using state.time in the g function

I am modelling an advection-diffusion-reaction model using the MATLAB PDE Toolbox. From reading the official documentation (Specify Boundary Conditions - MATLAB & Simulink - MathWorks Australia), the variable state.time is a scalar. When I use a static number for time (t=0.4) in the for loop, I don't have errors in my code but when I use state.time, I receive the following errors:

"

Unable to perform assignment because the left and right sides have a different number of elements.

Error in diffusion_equation>mygfun (line 249)

b(i) = y(tspan==state.time);

Error in pde.internal.pde2DBCSfnImpl/callNeumannFuncOnEdge (line 887)

bci = func(appRegion, state);

Error in pde.internal.pde2DBCSfnImpl/setNeumannBCOnEdge (line 576)

edgeG = self.callNeumannFuncOnEdge(bci,xyAllEdgeNodes, sPts, bci.g, ...

...

"

I want to include the time dependency of the bfun in gcoeff and I don't know how adjust my code to do so.

Thanks in advance.

function gcoeff = mygfun(~,state)
    Da = 0.01;
    
    [u_n, u_m] = size(state.u);
    b = zeros(u_n, u_m);
    tspan = 0:0.05:1;
    
    for i = 1:u_m 
        % This is solving for one governing equation, hence one         
        % for loop and the notation u(i). 
        bfun = @(t,y) state.u(i).*(1-y) - y./2;
        [t,y] = ode45(bfun, tspan, 0);

        b(i) = y(t==state.time);
%         b(i) =  y(t==0.4);
    end

    gcoeff = Da.*state.u.*(1-b)- b./2;

end

Side note: the vector tspan is the same vector used in the PDE solver.

Answer 1

[Zuber Khan]

Hi,

I understand that you are facing dimensionality error in your code.

It is possible that there is a mismatch in the dimensions of the variables you are working with in your mygfun function, specifically with the line

b(i) = y(tspan==state.time);

This suggests that y(tspan==state.time) is not returning a single element, which is causing the assignment error.

In order to debug and fix this issue, ensure that y(tspan==state.time) returns a single value. According to me, you can follow the approach as mentioned below to troubleshoot and resolve the problem:

1. Ensure that 'state.time' is a single scalar value and that 'tspan' contains the time points you expect.
2. Print or check the values of 'tspan' and 'state.time' to ensure they match up correctly.
3. If y(tspan==state.time) returns multiple values, you might need to adjust your logic to ensure a single value is assigned.

You can refer to the following modified mygfun function code snippet to debug and potentially fix the issue.

function gcoeff = mygfun(~, state)
    
    global y tspan
    Da = 0.01;
    
    [u_n, u_m] = size(state.u);
    b = zeros(u_n, u_m);
    tspan = 0:0.05:1;
    
   for i = 1:u_m 
        % This is solving for one governing equation, hence one         
        % for loop and the notation u(i). 
        bfun = @(t,y) state.u(i).*(1-y) - y./2;
        [t,y] = ode45(bfun, tspan, 0);
        
        % Print state.time and check against tspan
        disp(['state.time: ', num2str(state.time)]);
        disp(['tspan: ', num2str(tspan)]);

        % Find the index of the closest time point in tspan to state.time
        [~, idx] = min(abs(tspan - state.time));

        % Ensure that the index is valid and assign the corresponding value from y
        if ~isempty(idx) && idx > 0 && idx <= length(y)
            b(i) = y(idx);
        else
            error('No valid time point found in tspan for the given state.time.');
        end
   
   end
   
   gcoeff = Da.*state.u.*(1-b)- b./2;

end
   

In the above code, wherever (tspan==state.time) does not hold true, the closest possible value will get assigned to b(i). This is necessary to avoid any floating-point precision issues because if 'state.time' and the values in 'tspan' are not exactly the same due to numerical precision, the comparison (tspan==state.time) may not return 'true' for any element, leading to an empty result.

Hope it helps!

Regards

Zuber