How to include partial derivative term of matrix in finite difference equation?

I'm solving finite difference equation by tridiagonal system along time loop. I want to include integrodifferential term in my equation. For that purpose, I solved integration in a numerical way using trapezoidal rule and I got output values like matrix form. Now I want to solve the partial derivative of matrix(integral output) w.r.to x and include it in my equation. To solve partial derivative, I'm using GRADIENT function. but it makes an error: "Unable to perform assignment because the left and right sides have a different number of elements".

This is my code:

ymax=20; m=80; dy=ymax/m; %y=dy:dy:ymax;      %'i'th row
xmax=1; n=20; dx=xmax/n;                     %'j'th column   
tmax=100; nt=500; dt=tmax/nt; t=0:dt:tmax;   
phi=45; gr=10^6;
UWALL=zeros(m,nt); UOLD=zeros(m,nt); UNEW=0;
VOLD=zeros(m,nt); TNEW=0; TOLD=TNEW*ones(m,nt); TWALL=ones(1,length(t));
A=zeros([1,m]); 
B=A;
C=A;
D=A;
T=TOLD; U=UOLD; 
for j=1:nt
        for i=1:m
            if j>i
                C(i)=dt*VOLD(i,j)/4*dy-dt/(2*dy^2);
            elseif i>j
                A(i)=-dt*VOLD(i,j)/4*dy-dt/(2*dy^2);
            elseif i==j
                B(i)=1+dt*UOLD(i,j)/2*dx+dt/dy^2;
            end
        end
    end
    for j=2:nt
        for i=2:m-1
            if i==2
                D(i)=-(dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1)/2*dx))+(dt*(TWALL(j)-2*TOLD(i,j)+TOLD(i-1,j)/(2*dy^2)))-(dt*VOLD(i,j)*(TOLD(i-1,j)-TWALL(j)+TOLD(i,j)/4*dy))-((-dt)*(VOLD(i,j)/4*dy))-(dt*(TWALL(j)/2*dy^2));
            elseif i==m-1
                D(i)=-(dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1)/2*dx))+(dt*(TOLD(i+1,j)-2*TOLD(i,j)+TNEW/(2*dy^2)))-(dt*VOLD(i,j)*(TNEW-TOLD(i+1,j)+TOLD(i,j)/4*dy))-(dt*(VOLD(i,j)/4*dy))-(dt*(TNEW/2*dy^2));
            else
                D(i)=-(dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1)/2*dx))+(dt*(TOLD(i+1,j)-2*TOLD(i,j)+TOLD(i-1,j)/2*dy^2))-(dt*VOLD(i,j)*(TOLD(i-1,j)-TOLD(i+1,j)+TOLD(i,j)/4*dy));
            end
        end
        T(:,j)=TriDiag(A,B,C,D);             %call tridiagonal function
        dt=0.2+dt;
        TOLD=T;
    end
%============CALCULATE INTEGRAL==========================%
F=@(y) T;
a=0; b=20;
h=(b-a)/m;
i=1:1:m-1;
S=F(a+i*h);
I=(h./2)*(F(a)+2*sum(S)+F(b));
%================CALCULATE PARTIALDERIVATIVE=============%
%[IX,IY] = gradient(I);  %how to solve this derivative

    for j=1:nt
        for i=1:m
            if j>i
                C(i)=dt*VOLD(i,j)/4*dy-dt/(2*dy^2);
            elseif i>j
                A(i)=-dt*VOLD(i,j)/4*dy-dt/(2*dy^2);
            elseif i==j
                B(i)=1+dt*UOLD(i,j)/2*dx+dt/dy^2;
            end
        end
    end
    for j=2:nt
        for i=2:m-1
            if i==2                                    %here.......................................................................................want to include that derivative here.
                D(i)=-(dt*UOLD(i,j)*(-U(i,j-1)+UOLD(i,j)-UOLD(i,j-1)/2*dx))+(dt*(UWALL(j)-2*UOLD(i,j)+UOLD(i-1,j)/2*dy^2))+(dt*(gr^(-1/4)*cos(phi)*[IX,IY]+sin(phi)*(1/2)*(TWALL(j)+TOLD(i,j))))-(dt*VOLD(i,j)*(UOLD(i-1,j)-UWALL(j)+UOLD(i,j)/4*dy))-((-dt)*(VOLD(i,j)/4*dy))-(dt*(UWALL(j)/2*dy^2));
            elseif i==m-1
                D(i)=-(dt*UOLD(i,j)*(-U(i,j-1)+UOLD(i,j)-UOLD(i,j-1)/2*dx))+(dt*(UOLD(i+1,j)-2*UOLD(i,j)+UNEW/2*dy^2))+(dt*gr^(-1/4)*cos(phi)*[IX,IY]+sin(phi)*(1/2)*(TNEW+TOLD(i,j)))-(dt*VOLD(i,j)*(UNEW-UOLD(i+1,j)+UOLD(i,j)/4*dy))-(dt*(VOLD(i,j)/4*dy))-(dt*(UNEW/2*dy^2));
            else
                D(i)=-(dt*UOLD(i,j)*(-U(i,j-1)+UOLD(i,j)-UOLD(i,j-1)/2*dx))+(dt*(UOLD(i+1,j)-2*UOLD(i,j)+UOLD(i-1,j)/2*dy^2))-(dt*VOLD(i,j)*(UOLD(i-1,j)-UOLD(i+1,j)+UOLD(i,j)/4*dy))+(dt*(gr^(-1/4)*cos(phi)*[IX,IY]+sin(phi)*(1/2)*(T(i,j)+TOLD(i,j))));
            end
        end
        U(:,j)=TriDiag(A,B,C,D);
        dt=0.2+dt;
        UOLD=U;
    end

    
 %tridiagonal function
 
function x = TriDiag(e,f,g,r)
%input
% e = subdiagonal vector
% f = diagonal vector
% g = superdiagonal vector
% r = right hand side vector
% output:
% x = solution vector


n=length(f);
for k = 2:n
    factor = e(k)/f(k-1);
    f(k) = f(k) - factor*g(k-1);
    r(k) = r(k) - factor*r(k-1);
end
x(n) = r(n)/f(n);
for k = n-1:-1:1
    x(k) = (r(k)-g(k)*x(k+1))/f(k);
end
for i = 1:length(x)
    fprintf('\n%d = %f\n', i, x(i));
end

Please give some suggestion for solving derivative of matrix and how to include it in my equation?. Thank you for advance.

[Torsten]

Every question like yours should start with the problem formulation in a mathematical form: equation(s) to be solved, initial conditions and boundary conditions.

I hope you understand that it will be too hard for us to extract this information from your code.

[Nathi]

@Torsten Are you saying that we have to solve the partial derivative numerically like I solved the integral?.

[Torsten]

As far as I understand, all what you do above is one big problem. And since I don't understand what you mean by

Now I want to solve the partial derivative of matrix(integral output) w.r.to x and include it in my equation.

and how it is part of the overall problem, it would help if you included the equations, boundary and initial conditions and explained which part of this problem is related to the task you try to get an answer for.

[Nathi]

@Torsten The issue of my code is momentum equation which has partial derivative term to run the overall code. Firstly, my code is based on the Governing equation and it contains continuity, momentum(velocity) i.e.U and energy(temperature) i.e.T along initial and boundary condition and I'm using this equations in a numerical solution, especially implicit crank nicolson scheme. " The process of this method is by solving energy equation only I can solve momentum equation and this equation has integrodifferential part (as of now the issue of the term) and then the continuity equation".

And I solved energy equation in the system of tridiagonal matrix since I'm using crank nicolson scheme and the output value is T which I got 2D matrix form

for j=1:nt                %energy equation
        for i=1:m
            if j>i
                C(i)=dt*VOLD(i,j)/4*dy-dt/(2*dy^2);
            elseif i>j
                A(i)=-dt*VOLD(i,j)/4*dy-dt/(2*dy^2);
            elseif i==j
                B(i)=1+dt*UOLD(i,j)/2*dx+dt/dy^2;
            end
        end
    end
    for j=2:nt
        for i=2:m-1
            if i==2
                D(i)=-(dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1)/2*dx))+(dt*(TWALL(j)-2*TOLD(i,j)+TOLD(i-1,j)/(2*dy^2)))-(dt*VOLD(i,j)*(TOLD(i-1,j)-TWALL(j)+TOLD(i,j)/4*dy))-((-dt)*(VOLD(i,j)/4*dy))-(dt*(TWALL(j)/2*dy^2));
            elseif i==m-1
                D(i)=-(dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1)/2*dx))+(dt*(TOLD(i+1,j)-2*TOLD(i,j)+TNEW/(2*dy^2)))-(dt*VOLD(i,j)*(TNEW-TOLD(i+1,j)+TOLD(i,j)/4*dy))-(dt*(VOLD(i,j)/4*dy))-(dt*(TNEW/2*dy^2));
            else
                D(i)=-(dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1)/2*dx))+(dt*(TOLD(i+1,j)-2*TOLD(i,j)+TOLD(i-1,j)/2*dy^2))-(dt*VOLD(i,j)*(TOLD(i-1,j)-TOLD(i+1,j)+TOLD(i,j)/4*dy));
            end
        end
        T(:,j)=TriDiag(A,B,C,D);             %call tridiagonal function
        dt=0.2+dt;
        TOLD=T;        %  got output values 'T' here
    end

and then now I have to solve momentum equation. however I couldn't to solve this equation because it has partial derivative with intergral term which depends on T.

For that purpose first i solved integral part using trapezoidal rule because i'm dealing with numerical solution.ie.,

  %============CALCULATE INTEGRAL==========================%
F=@(y) T;                       % I took T which is the output value of energy as a function w.r.to Y 
a=0; b=20;
h=(b-a)/m;
i=1:1:m-1;
S=F(a+i*h);
I=(h./2)*(F(a)+2*sum(S)+F(b));  

and then I want to solve partial derivative of( I=(h./2)*(F(a)+2*sum(S)+F(b)); ) w.r.to x using Numerically.

that why I'm using gradient funtion

%================CALCULATE PARTIALDERIVATIVE=============%
[IX,IY] = gradient(I);  %how to solve this derivative and if i include this i got error. "unable to perform.....".

and next this [IX,IY] will include in momentum equation

    for j=1:nt
        for i=1:m
            if j>i
                C(i)=dt*VOLD(i,j)/4*dy-dt/(2*dy^2);
            elseif i>j
                A(i)=-dt*VOLD(i,j)/4*dy-dt/(2*dy^2);
            elseif i==j
                B(i)=1+dt*UOLD(i,j)/2*dx+dt/dy^2;
            end
        end
    end
    for j=2:nt
        for i=2:m-1
            if i==2            %here.......................................want to include that derivative here.
                D(i)=-(dt*UOLD(i,j)*(-U(i,j-1)+UOLD(i,j)-UOLD(i,j-1)/2*dx))+(dt*(UWALL(j)-2*UOLD(i,j)+UOLD(i-1,j)/2*dy^2))+(dt*(gr^(-1/4)*cos(phi)*[IX,IY]+sin(phi)*(1/2)*(TWALL(j)+TOLD(i,j))))-(dt*VOLD(i,j)*(UOLD(i-1,j)-UWALL(j)+UOLD(i,j)/4*dy))-((-dt)*(VOLD(i,j)/4*dy))-(dt*(UWALL(j)/2*dy^2));
            elseif i==m-1
                D(i)=-(dt*UOLD(i,j)*(-U(i,j-1)+UOLD(i,j)-UOLD(i,j-1)/2*dx))+(dt*(UOLD(i+1,j)-2*UOLD(i,j)+UNEW/2*dy^2))+(dt*gr^(-1/4)*cos(phi)*[IX,IY]+sin(phi)*(1/2)*(TNEW+TOLD(i,j)))-(dt*VOLD(i,j)*(UNEW-UOLD(i+1,j)+UOLD(i,j)/4*dy))-(dt*(VOLD(i,j)/4*dy))-(dt*(UNEW/2*dy^2));
            else
                D(i)=-(dt*UOLD(i,j)*(-U(i,j-1)+UOLD(i,j)-UOLD(i,j-1)/2*dx))+(dt*(UOLD(i+1,j)-2*UOLD(i,j)+UOLD(i-1,j)/2*dy^2))-(dt*VOLD(i,j)*(UOLD(i-1,j)-UOLD(i+1,j)+UOLD(i,j)/4*dy))+(dt*(gr^(-1/4)*cos(phi)*[IX,IY]+sin(phi)*(1/2)*(T(i,j)+TOLD(i,j))));
            end
        end
        U(:,j)=TriDiag(A,B,C,D);
        dt=0.2+dt;
        UOLD=U;
    end

The major line of the momentum equation which i want to include is

    if i==2            %here......................................................................................want to include that derivative here.
                D(i)=-(dt*UOLD(i,j)*(-U(i,j-1)+UOLD(i,j)-UOLD(i,j-1)/2*dx))+(dt*(UWALL(j)-2*UOLD(i,j)+UOLD(i-1,j)/2*dy^2))+(dt*(gr^(-1/4)*cos(phi)*[IX,IY]+sin(phi)*(1/2)*(TWALL(j)+TOLD(i,j))))-(dt*VOLD(i,j)*(UOLD(i-1,j)-UWALL(j)+UOLD(i,j)/4*dy))-((-dt)*(VOLD(i,j)/4*dy))-(dt*(UWALL(j)/2*dy^2));

This line not accepting the derivatie form.

I hope I explained it clearly now. and please try to understand it and clear this issue. This diminutive issue affected my overall code. once again thank you.

[Torsten]

I don't understand your equations in the way you present them. Sorry.

[Nathi]

@Torsten Still why you are not understand my equations?. I explained it very clearly.

once again see my full equation.

I splitting my equation for tridiagonal system and using A,B,C,D in my code. Here A,B,C are superdiagonal, diagonal and subdiagonal respectively and D is the Rhs of tridiagonal matrix. so I'm taking three multiple D equation such as upper wall, lower wall along boundary condition and interior nodes.

Please try it again. its very small issue once you understand my code or tell me where you didn't get it.

[Torsten]

And where do you solve the energy equation for temperature T ? Or what is T and how is it computed ?

Same for V. Is it given or is there an extra equation to be solved for it ?

[Nathi]

I solved the energy equation at the beginning for temperature T. The input value of T = TOLD And then the output value stored in T only which I called tridiagonal form like T(:,j)=TriDiag(A,B,C,D); (see the below code at the end line.)

%==sapce============
ymax=20; m=80; dy=ymax/m; %y=dy:dy:ymax;      %'i'th row
xmax=1; n=20; dx=xmax/n;                     %'j'th column
%=time========
tmax=100; nt=500; dt=tmax/nt; t=0:dt:tmax; 

phi=45; gr=10^6;
%============initialization=============
UWALL=zeros(m,nt); UOLD=zeros(m,nt); UNEW=0;
VOLD=zeros(m,nt); TNEW=0; TOLD=TNEW*ones(m,nt); TWALL=ones(1,length(t));
A=zeros([1,m]); 
B=A;
C=A;
D=A;
T=TOLD; U=UOLD; 
%=======ENERGY===========================
for j=1:nt
        for i=1:m
            if j>i
                C(i)=dt*VOLD(i,j)/4*dy-dt/(2*dy^2);
            elseif i>j
                A(i)=-dt*VOLD(i,j)/4*dy-dt/(2*dy^2);
            elseif i==j
                B(i)=1+dt*UOLD(i,j)/2*dx+dt/dy^2;
            end
        end
    end
    for j=2:nt
        for i=2:m-1
            if i==2
                D(i)=-(dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1)/2*dx))+(dt*(TWALL(j)-2*TOLD(i,j)+TOLD(i-1,j)/(2*dy^2)))-(dt*VOLD(i,j)*(TOLD(i-1,j)-TWALL(j)+TOLD(i,j)/4*dy))-((-dt)*(VOLD(i,j)/4*dy))-(dt*(TWALL(j)/2*dy^2));
            elseif i==m-1
                D(i)=-(dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1)/2*dx))+(dt*(TOLD(i+1,j)-2*TOLD(i,j)+TNEW/(2*dy^2)))-(dt*VOLD(i,j)*(TNEW-TOLD(i+1,j)+TOLD(i,j)/4*dy))-(dt*(VOLD(i,j)/4*dy))-(dt*(TNEW/2*dy^2));
            else
                D(i)=-(dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1)/2*dx))+(dt*(TOLD(i+1,j)-2*TOLD(i,j)+TOLD(i-1,j)/2*dy^2))-(dt*VOLD(i,j)*(TOLD(i-1,j)-TOLD(i+1,j)+TOLD(i,j)/4*dy));
            end
        end
        T(:,j)=TriDiag(A,B,C,D);             %call tridiagonal function
        dt=0.2+dt;
        TOLD=T;
    end
..........velocity equation......................

In addition, There is extra equation to be solved i.e. continuity equation for V. I didn't include the continuity equation yet. Because before solving 'V' I have to solve momentum equation 'U'.

[Torsten]

It's not clear from your notation, but my guess is that for each value of j, you have to evaluate

integral_{y(j)}^{Inf}  T(x(i),y) dy

which can be approximated as

I(i,j) = trapz(y(j:end),T(i,j:end))

For j fixed, this gives you an x-profile I(1,j), I(2,j),...,I(end,j).

Now apply gradient on this vector:

gradient(I(:,j))/gradient(x(:))

Maybe one advice (although I know you most probably will not follow it): First make the pure flow equations work, then you can try to include extra terms from the energy equation.

[Nathi]

Yes! I made it. This is the baseform only. For the system, I need to solve that integral term in momentum equation. That's why I try to solve this case.