diffrential equation using ode45

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

[Walter Roberson]

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]

dont take xm

[abhishek singh]

please get it output in matrix form [x(:1),x(:2);x(:3),x(:4)]

[abhishek singh]

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]

I do not understand what you mean by "dont take xm" ?

Please show your ode45() call.

[abhishek singh]

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]

If you are using xm then w=1+xm

[abhishek singh]

Sir you understand the problem

[Walter Roberson]

What is ou_seq ?

[abhishek singh]

this is random noise i am using this in w =1+xm. xm is functional form of noise

[Walter Roberson]

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]

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).

Answer 1

[Walter Roberson]

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))];

[abhishek singh]

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]

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]

trace of any density matrix is always one .

[abhishek singh]

Is there any other way to add noise without using interpolation method in this code .

[abhishek singh]

The code you send is not properly working.

[Walter Roberson]

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]

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.