Family of ODE with summation

I need to solve this family of ODE with summation inside and I think I'm not using properly ODE45. I really need some help.

Where K_Br is a matrix of known values previously calculated.

I wrote a function to define the system of ODE:

function ndot = System_ni (t,n)
ndot=zeros(M,1);
n=zeros(M,1);
V=zeros(M,1);

    for i = 1:M
        for j=1:i-1
            V(j)=K_Br(j,i-j).*n(j);
        end
        nw(i)=1/2*sum(V(1:i-1).*n(1:i-1));
        for j=1:M
            VV(j)=K_Br(i,j).*n(j);
        end
        nnw(i)=n(i)*sum(VV(1:M));
        ndot(i) =nw(i)+nnw(i);
    end
    ndot(1)=n(1)*sum(K_Br(1,:)*n(:));
end

And I'm recalling in in ODE45 with the intial condition.

cond_in=[N_in;zeros(M-1,1)];
[t,n] = ode45(@System_ni,[0 T],cond_in);

But something it's not working and I think I'm making some mistakes with ODE45.

Thank you in advance for any help!

[Walter Roberson]

It is not obvious to me why you are overwriting the already-calculated ndot(1) at the end ?

[Sara Crosetto]

Because equation is different for ndot(1) since first summation is zero.

Actually, I should have written it at the very beginning, I think.

Answer 1

[Alan Stevens]

Your System_ni function doesn't seem to have acces to your K_Br matrix. You need to pass it into System_ni or specify it within System_ni.

[Sara Crosetto]

I've specified it inside System_ni, I didn't report it for the sake of brevity. The whole function actually is:

function ndot = System_ni (t,n)
Lo=1.68e-07;
M=500;
L=20*Lo;
l=linspace(Lo,L,M);

Kb=1.38e-23;
eta_w=9e-04;
Temp=337.15;
C=(2*Kb*Temp)/(3*eta_w);

K_Br=zeros(M,M);
for i=1:M
    for j=1:M
    K_Br(i,j)=C*((l(i)+l(j))^2/(l(i)*l(j)));
    end 
end

ndot=zeros(M,1);
n=zeros(M,1);
V=zeros(M,1);

    for i = 1:M
        for j=1:i-1
            V(j)=K_Br(j,i-j).*n(j);
        end
        nw(i)=1/2*sum(V(1:i-1).*n(1:i-1));
        for j=1:M
            VV(j)=K_Br(i,j).*n(j);
        end
        nnw(i)=n(i)*sum(VV(1:M));
        ndot(i) =nw(i)+nnw(i);
    end
    ndot(1)=n(1)*sum(K_Br(1,:)*n(:));
end

[Alan Stevens]

You pass in n, but then overwrite it with zeros

n=zeros(M,1);

[Sara Crosetto]

I removed it, but it doesn't work anyway..

[Alan Stevens]

What is the error message? If you upload all your code, we could investigate in more detail.

[Sara Crosetto]

The error is:

Warning: Failure at t=1.066340e+01.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.842171e-14) at
time t. 

Solution is not correct for time <t=1.066340e+01, anyway.

The whole code and the function are attached. Thank you so much for your kind help!

[Alan Stevens]

Try the following (note: I've included the System_ni function at the end of the Report8 script. You will need to decide if the end result is sensible.

N_in=1.7521e+14;

N=100;
T=500;
Lo=1.68e-07;
M=500;
L=20*Lo;
l=linspace(Lo,L,M);

Kb=1.38e-23;
eta_w=9e-04;
Temp=337.15;
C=(2*Kb*Temp)/(3*eta_w);

K_Br=zeros(M,M);
for i=1:M
    for j=1:M
    K_Br(i,j)=C*((l(i)+l(j))^2/(l(i)*l(j)));
    end 
end

cond_in=[N_in;zeros(M-1,1)];
[t,n] = ode45(@(t,n) System_ni(t,n,K_Br),[0 T],cond_in);  % Use this form to pass an extra parameter to System_ni
plot(t,n(:,1))


function ndot = System_ni (~,n,K_Br)  % Take K_Br from previous calcuation rather than recalculating it.
M = size(K_Br,1);
ndot=zeros(M,1);
V=zeros(M,1);
nnw=zeros(M,1);
nw=zeros(M,1);

    for i = 1:M
        for j=1:i-1
            V(j)=K_Br(j,i-j).*n(j);
        end
        nw(i)=1/2*sum(V(1:i-1).*n(1:i-1));
        VV = K_Br(i,:)*n;           % row vector * column vector automatically does summation
        nnw(i)=n(i)*VV;             % VV is a single value
        ndot(i) = nw(i) - nnw(i);   % nnw needs negative sign
    end
    
end