diffrential equation using ode45

2019 12 18
1 个回答
更新时间：2019 12 21
10 次查看（30 天）

0 个投票

how to solve this diffrential equation using ode45
function dxdt = decoherence1nv(t,x,del,delt)
xm = interp1(delt,del,t);
dxdt = zeros(4,1);
a=0.04; Q=0.5;e=0.1;z=pi/3;w=1;
dxdt(1) = -a .* x(1) + 1i .*0.5.* e .* Q .* sin(z) .* x(2) - 1i.*0.5.* e .* Q .* sin(z) .* x(3) ;
dxdt(2) = 1i .*0.5.* e .* Q .* sin(z).* x(1) - 1i .* w .* x(2) - 0.5 .* a .* x(2) -1i .* e .* Q .* cos(z) .* x(2) - 1i .*0.5.* e .* Q .* sin(z) .* x(4) ;
dxdt(3) = -1i .*0.5* e .* Q .* sin(z) .* x(1) + 1i .* w .* x(3) - 0.5 .* a .* x(3) + 1i .* e .* Q .* cos(z) .* x(3) + 1i.*0.5* e .* Q .* sin(z) .* x(4) ;
dxdt(4) = a .* x(1) - 1i .*0.5* e .* Q .* sin(z) .* x(2) + 1i .*0.5* e .* Q .* sin(z) .* x(3);
end

12 个评论
显示 10更早的评论 隐藏 10更早的评论

Walter Roberson
Walter Roberson 2019-12-18
What problem are you running into?
xm = interp1(delt,del,t);
That would typically introduce discontinuities in the calculations. You should be breaking the single ode45() call over the timespan, into one call for each different boundary of delt, so that you would not have discontinuities in the second derivative.
Or rather you would do that if not for the fact that you do not use xm in your calculation.
abhishek singh
abhishek singh 2019-12-18
please get it output in matrix form [x(:1),x(:2);x(:3),x(:4)]
abhishek singh
abhishek singh 2019-12-18
this i am taking as xm which is used in w=1+xm
function x=ou_seq(nt,tmax,tau,c)
%_start=0;%simulation start time
%t_end=10;%simulation end time
%dt=0.01;
%tau=.05;%relaxation
%c=1;%diffusion constant
x0=0; %initial value of stochastic variable
%mu=0;%mean of stochastic process x
%y0=0; %initial value for integral x
%start_dist=-2.0; %start of OU pdf
%end_dist=2.0; %end of OU pdf
%trange=[0 tmax];
tmax = 100;
nt=100;tau=.5; c=10;
T=linspace(0,tmax,nt);
dt=tmax/(nt-1);
%variance=(c*tau*0.5)*(1-(exp(-dt/tau))^2);
%compute x and y
i=1;
x(1)=x0;
%y(1)=y0;
for t=T(1)+dt:dt:T(nt)
i=i+1;
r1=randn;
% r2=randn;
x(i)=x(i-1)*exp(-dt/tau)+sqrt((c*tau*0.5)*(1-(exp(-dt/tau))^2))*r1;
end
%plot(T,x)
Walter Roberson
Walter Roberson 2019-12-18
I do not understand what you mean by "dont take xm" ?
Please show your ode45() call.
abhishek singh
abhishek singh 2019-12-18
clear all
% Define L and zrange
% this program has been tested on Feb. 3, 2014 for its compatibility with
%earlier program with fixed Delta. This works fine for given fluctuating
%random field supplied bi xia.
%Delta + or - does not matter. Only net magnitude of Delta matters here.
tic;
format long
tmax = 10;
nt=100;
itrm=1000;
%eps=1e-10;
tau=.5; c=10;
dt=tmax/(nt-1);
x1 = 1;x2 = 0;x3 = 0;x4 = 0;% set initial condition
x0=[x1 x2 x3 x4]; % intial conditions matrix
delt=linspace(0,tmax,nt);
for itr=1:itrm
xia=ou_seq(nt,tmax,tau,c);
%xia=0*delt;
del=xia;
[t,x] = ode113(@(t,x)decoherence1nv(t,x,del,delt),0:dt:tmax,x0);
%[t,x] = ode113(@(t,x)decoherence1nv(t,x),0:dt:tmax,x0);
rho11=x(:,1);
rho12=x(:,2);
rho21=conj(x(:,2));
rho22=x(:,4);
%puret=zeros(itrm,nt);
%entropy=zeros(itrm,nt);
for it=1:length(t)
rhot=[rho11(it),rho12(it);rho21(it),rho22(it)];
sqrhoc=rhot*rhot;
ea=eig(rhot);
%sum(ea)
%entropy(itr,it)=sum(-ea.*(log(ea+eps)/log(2)));
puret(itr,it)=trace(sqrhoc);
%entropy(it)=sum(-ea.*(log(ea+eps)/log(2)));
%puret(it)=trace(sqrhoc);
%clear ea rhot sqrhoc
end
trace(rhot)
end
%trace(rhot)
puretf=mean(puret,1);
%entrof=mean(entropy,1);
hold on
figure(1)
plot(t,puretf,'r')
%figure(2)
%plot(t,entrof,'k')
%save paper_fig_constant_field.mat
%%%case3 is used in general
%save nonoise_itr_10_tx_10_2nd_case5.mat
%save constant_driving_itr_1000_tmax_10_case1.mat
toc;
abhishek singh
abhishek singh 2019-12-18
If you are using xm then w=1+xm
abhishek singh
abhishek singh 2019-12-18
Sir you understand the problem
abhishek singh
abhishek singh 2019-12-18
this is random noise i am using this in w =1+xm. xm is functional form of noise
Walter Roberson
Walter Roberson 2019-12-18
I do not understand the problem. The code does not produce any error messages when it is run, and after the code, x is already in the form [x(:1),x(:2);x(:3),x(:4)] so I do not know what you are looking for.
Your ou_seq function should not be ignoring its inputs.
Using interp1() to calculate xm would be a discontinuity problem in your decoherence1nv function, except that that function ignores xm after it calculates it.
The ode*() series of function calls estimate gradient and hessian, and in order to do that, the function behaviour must be continuous to two derivatives within the time range that is being calculated over. If you use w=1+xm then you would be violating that condition at every delt boundary. You need to rewrite the code to make separate ode*() calls over time range delt(1:2), delt(2:3), delt(3:4) and so on.
abhishek singh
abhishek singh 2019-12-18
how can i write the seprate code for (You need to rewrite the code to make separate ode*() calls over time range delt(1:2), delt(2:3), delt(3:4) and so on).

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2019-12-18

0 个投票

current_x0 = x0;
for idx = 1:length(delt)-1
ts = delt(idx); te = delt(idx+1); tm = (ts+te)/2;
[part_t{idx},part_x{idx}] = ode113(@(t,x)decoherence1nv(t,x,del,delt),[ts,dm,te],current_x0);
current_x0 = part_x{idx}(end,:);
end
SkipFirst = @(V) V(2:end,:);
t = [part_t{1}{1}; cell2mat(cellfun(SkipFirst, part_t, 'uniform', 0))];
x = [part_x{1}(1,:); cell2mat(cellfun(SkipFirst, part_x, 'uniform', 0))];

7 个评论
显示 5更早的评论 隐藏 5更早的评论

abhishek singh
abhishek singh 2019-12-18
sir ,accutually for a=0, trace(rhot) must be equal to one (trace(rhot)=1) and trace(rhot*rhot) also equal to one but it exceed from one i dont understand why?
[t,x] = ode113(@(t,x)decoherence1nv(t,x),0:dt:tmax,x0);
rho11=x(:,1);
rho12=x(:,2);
rho21=x(:,3);
rho22=x(:,4);
%puret=zeros(itrm,nt);
%entropy=zeros(itrm,nt);
for it=1:length(t)
rhot=[rho11(it),rho12(it);rho21(it),rho22(it)];
sqrhoc=rhot*rhot;
plot(t,sqrhoc)
function dxdt = decoherence1nv(t,x)
dxdt = zeros(4,1);
a=0.000; Q=0.5;e=0.1;z=pi/3;w=1+xm;
dxdt(1) = -a .* x(1) + 1i .*0.5.* e .* Q .* sin(z) .* x(2) - 1i.*0.5.* e .* Q .* sin(z) .* x(3) ;
dxdt(2) = 1i .*0.5.* e .* Q .* sin(z).* x(1) - 1i .* w .* x(2) - 0.5 .* a .* x(2) -1i .* e .* Q .* cos(z) .* x(2) - 1i .*0.5.* e .* Q .* sin(z) .* x(4) ;
dxdt(3) = -1i .*0.5* e .* Q .* sin(z) .* x(1) + 1i .* w .* x(3) - 0.5 .* a .* x(3) + 1i .* e .* Q .* cos(z) .* x(3) + 1i.*0.5* e .* Q .* sin(z) .* x(4) ;
dxdt(4) = a .* x(1) - 1i .*0.5* e .* Q .* sin(z) .* x(2) + 1i .*0.5* e .* Q .* sin(z) .* x(3);
end
Walter Roberson
Walter Roberson 2019-12-18
ode113() is permitted to produce invalid results when you ask it to calculate across a discontinuity.
At the moment, I do not see any mathematical reason why trace(rht) should ever be 1.
abhishek singh
abhishek singh 2019-12-18
trace of any density matrix is always one .
abhishek singh
abhishek singh 2019-12-21
Is there any other way to add noise without using interpolation method in this code .
abhishek singh
abhishek singh 2019-12-21
The code you send is not properly working.
Walter Roberson
Walter Roberson 2019-12-21
Is there any other way to add noise without using interpolation method in this code .
Yes, there is.
As I indicated before, the interpolation you are doing introduces discontinuities in the Hessian, and because that is not permitted, you need to break the calculation up into segments in which there is no discontinuity.
With there not being any discontinuity, instead of using interp1(), you can pass in the information you need to do the interpolation yourself. You would need to know the base time of this segment, and the base noise, xm(K), and the slope m=(xm(K+1)-xm(K))/delta_t . The interpolation would then be (t-base_t) * slope + base_noise, which would be faster than calling interp1()
Walter Roberson
Walter Roberson 2019-12-21
Note also that in theory you could fit a degree 4+2=6 or higher polynomial to the noise and use the fitted data with interp1(): if you were to do that, then the theoretical problems about the noise giving discontinuous hessians would not apply. However, a degree N polynomial would of course have at most ceil(N/2) peaks, and degree 7 is probably about the maximum polynomial order that is reasonably numerically stable, so I doubt that it would be an acceptable model of your physics to fit your noise to degree 6 to 8, so you will probably need to split the timespan like I discussed earlier.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Correlation and Convolution 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by