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
2 次查看(过去 30 天)
显示 更早的评论
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
2 个评论
David Goodmanson
2018-1-6
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?
采纳的回答
David Goodmanson
2018-1-7
编辑:David Goodmanson
2018-1-7
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
8 个评论
David Goodmanson
2018-1-18
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!