Problem solving a differential equation with time dependent parameter.
1 次查看(过去 30 天)
显示 更早的评论
I am trying to solve a system of rate equations, with time dependent ionization rate w(t) (my time dependent parameters in my case): my time vector is -10:0.1:10 so w(t) is also a vector and my equations are
dN0/dt=-W(t)*N0(t) and dN1/dt=W(t)*N0(t)
Think of it as a system with two levels, and the population at each level for time t is N0(t)and N1(t). An alternative approach for the equation would be to use the dN/dt=w(t)[Ng-N(t)] where Ng is a number, the initial population in ground state.
I am trying to solve this problem using the ode functions. However, using ode45 is not possible, since it takes a lot of time to compute and I cannot compute it and the other solvers don't give the "expected" solution.
I define my function as:
function dNdt=rate_equations(t,y,ft,f)
f = interp1(ft,f,t);
% dNdt=[-f.*y(1) ; f.*y(1)];
dNdt=f.*[2.4*1e19-y];
end
the comment denotes another alternative approach and in the main script I calculate the solution as:
initial_cond=[0]; %initial_cond=[2.4*1e19;0]; when using the system pair
t_span=-10:0.1:10;
a=wt';
b=W_time';
[tt,NN]=ode23t(@(t,y)rate_equations(t, y, a, b), t_span, initial_cond);
My time dependent parameter vector W_time is defined as:
E_a=5.14*1e9 ;
omega_au=4.134*1e16;
eta0=376.6;
E_omega=sqrt(2*eta0*1e14); %Field of fundamental
E_omega2=sqrt(2*eta0*2e13);% field of second harmonic
theta=pi/2;
wt=-10:0.01:10 ;
W_time=4.*omega_au.*(E_a./abs(E_omega.*cos(wt)+E_omega2.*cos(2.*wt+theta))).*exp((-2/3).*(E_a./abs(E_omega.*cos(wt)+E_omega2.*cos(2.*wt+theta))))
I don't get the results I expected. In theory I should be getting something like this:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/145754/image.jpeg)
but with the plot saturating at some point. When I use the stiff solvers I get some very slight oscillations. The non stiff ones and ode45 are calculating for ever.
So my question is, what I am doing wrong? And how can I improve the computing time of the ode45 solver?
回答(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!