hi everyone,I write a code for solving N differential first order coupled equations ,I want to know if there is better way to do it , or more efficient way , thanks

the equations:

(dC_k)/dt= -i∑C_n (t) *R_nk *(e^(-i(ω_n-ω_k-ω)t)+e^(-i(ω_n-ω_k+ω)t))

%the sum is over n 

K=1:N

n=1:N

R  is a matrix of numbers , R_nk is R(n,k) ,for different n there is different  ω_n  ,ω is  
constant frequency ,ω_n,ω_k  are  frequencies

my code:

[t,c] = ode45(@(t,c)myode(t,c,cof,energy,N,w),[0 tf],initial_cond);

function dc = myode(t,c,cof,energy,N,w)

%N coupled differential equations 

phase_minus=@(w_nk,w) exp(-1i*(w_nk-w));

phase_plus=@(w_nk,w)  exp(-1i*(w_nk+w));

dc=zeros(N,1);

for k=1:N

    w_k=energy(k);

       for n=1:N   

       w_n=energy(n);

       w_nk=w_n-w_k;

       R_n=cof(k,n);

       dc(k)=dc(k)-1i*c(n)*(phase_plus(w_nk,w)^t+phase_minus(w_nk,w)^t)*R_n;

       end

end

end

[David Goodmanson]

Hi fatema,

Whether or not a different way works better may depend on the number of components of c. What is a typical value for N? Also, the code implies that R depends just on one index n and is not a matrix R(k,n), is that correct?

[fatema hamodi]

hi .it's my fault I should write R_nk rather than R_n .the number of the equations is large say 1000 , that's why I am looking for efficient way .

Answer 1

[David Goodmanson]

Hi fatema,

Here is an alternate ode function. It assumes that cof is NxN, c is a column vector of length N, and energy is a vector of length N. The code factors out the cos(w*t) part, sets up a matrix of values (w_n-w_k) (a difference matrix for the energies), does element-by-element multiplication with cof and then usual matrix multiplication times c.

For N = 1000, the run time goes from about 4 sec to .05 sec.

function dc = myode1(t,c,cof,energy,N,w)

energy = energy(:);                  % make sure it's a column vector
ww = -energy + energy.';
dc = (-2i)*cos(w*t)*(cof.*(exp(-i*ww*t)))*c;

% if you do not have a newer version of Matlab with implicit expansion, use
%  ww = -repmat(energy,1,N) + repmat(energy.',N,1);

end

[fatema hamodi]

 David Goodmanson thanks a lot  !!

[fatema hamodi]

 David Goodmanson my code and the code you suggested gives different numerical results , but I think mine is right but I am not sure , do you have any idea why the codes don't give the same results

[David Goodmanson]

Hi fatema,

I generally check the new version vs. the original version, and in this case I compared outputs of the two myode functions, figuring that if those two agreed then ode45 would give the same output. Here is an example:

% make up some data
N = 1000;
cof = rand(N,N)+i*rand(N,N);
c = rand(N,1)+i*rand(N,1);
energy = rand(1,N);
t = 1;
w = 1/2;

% try it out.  
dc  =  myode(t,c,cof,energy,N,w);    % original function
dc1 = myode1(t,c,cof,energy,N,w);    % new function from the answer

% find largest relative error of any component of dc1 compared to dc
max(abs((dc1-dc)./dc))
ans = 1.4492e-15

If you have data that does not work like this, I would be interested to see it.

[fatema hamodi]

hi I really appreciate your help thanks , I try your code and I was

astonished that they give the same result , in my system they don't .

try this code and you will see the difference for other systems they

are absolutely different , no matter the number of equations N (small or

large)

N = 2;

cof = [0.5449,0.866;0.866,0.0551];

energy = [-10.899;-1.1010];

w = 9.798;

tf=100;

initial_cond=zeros(1,N);

initial_cond(1)=1;

% try it out.

[t_1,y_1] = ode45(@(t,y)myode_old(t,y,cof,energy,N,w),[0 tf],initial_cond);

[t,y] = ode45(@(t,y)myode11(t,y,cof,energy,w),[0 tf],initial_cond);

figure

plot(t_1,abs(y_1(:,1).^2))

ylabel('numerical sol 1')

figure

plot(t,abs (y(:,1)).^2)

ylabel('numerical sol 2')

[fatema hamodi]

Furthermore I am not sure which method is right , that's why I really interested to know why the results are different !

[fatema hamodi]

I found out the reason, and it seems so bizarre for me , in matlab there is difference between cos(w*t) and 0.5*(exp(iwt)+exp(-iwt)) once you replace the two exponents in the original code , they will gives the same result , do you have any idea why is that ?

[David Goodmanson]

Hi fatema,

Since you are citing a mathematical identity, I don't think that's actually the reason. One of the things I noticed about your code before was that you have

phase_minus=@(w_nk,w) exp(-1i*(w_nk-w));
phase_plus=@(w_nk,w)  exp(-1i*(w_nk+w));
and
dc(k)=dc(k)-i*c(n)*(phase_plus(w_nk,w)^t+phase_minus(w_nk,w)^t)*R_n;

so you are using (exp(i*w)) ^ t. That is not a very good way to be doing things, and it is cleaner, better and faster to use exp(i*w*t) instead. Algebraically they should be the same so I didn't comment on that at the time, but they give different results somehow in ode45.

You could replace the original code with

phase_minus = @(w_nk,w) exp(-1i*(w_nk-w)*t);
phase_plus  = @(w_nk,w) exp(-1i*(w_nk+w)*t);
and
dc(k)=dc(k)-1i*c(n)*(phase_plus(w_nk,w)+phase_minus(w_nk,w))*R_n;

in which case the modified plot makes better sense than the old one, and also agrees with the myode1 plot.

The inline functions are short enough that I like not bothering with them and just going with

   dc(k)=dc(k)-1i*c(n)*(exp(-i*(w_nk -w)*t) + exp(-i*(w_nk +w)*t))*R_n;

In any case myode1 will be faster for large N. Let me know if you still do not get agreement in the methods.

[fatema hamodi]

now I get an agreement ,this help me so much ! , thanks a lot !