Error in ode45 solution

Ode45 results dont maintain power conservation for some cases. How to ensure that power conservation is maintained?

While using ode45 to solve coupled differential equation, when one of the input parameters is real,the solution to ode45 maintains power conservation. When i change that input parameter from real to complex then ode45 solution is not maintaing power conservation. I have used relatve error tolerance = 1e-6 and absolute error tolerance = 1e-9. Please let me know what can be done to get the correct results.

For eg, In the below code, tv* parameter was initially real, which maintained power conservation, however when i change tv* to complex, the solution doesnt maintain power conservation.

options=odeset('RelTol',1e-6,'AbsTol',1e-9);
f=@(z,xx_val) -1i*[tv12*xx_val(2)*exp(1i*del_beta12*z)+tv13*xx_val(3)*exp(1i*del_beta13*z)+tv14*xx_val(4)*exp(1i*del_beta14*z)+tv15*xx_val(5)*exp(1i*del_beta15*z)+tv16*xx_val(6)*exp(1i*del_beta16*z);...
    tv21*xx_val(1)*exp(1i*del_beta21*z)+tv23*xx_val(3)*exp(1i*del_beta23*z)+tv24*xx_val(4)*exp(1i*del_beta24*z)+tv25*xx_val(5)*exp(1i*del_beta25*z)+tv26*xx_val(6)*exp(1i*del_beta26*z);...
    tv31*xx_val(1)*exp(1i*del_beta31*z)+tv32*xx_val(2)*exp(1i*del_beta32*z)+tv34*xx_val(4)*exp(1i*del_beta34*z)+tv35*xx_val(5)*exp(1i*del_beta35*z)+tv36*xx_val(6)*exp(1i*del_beta36*z);...
    tv41*xx_val(1)*exp(1i*del_beta41*z)+tv42*xx_val(2)*exp(1i*del_beta42*z)+tv43*xx_val(3)*exp(1i*del_beta43*z)+tv45*xx_val(5)*exp(1i*del_beta45*z)+tv46*xx_val(6)*exp(1i*del_beta46*z);...
    tv51*xx_val(1)*exp(1i*del_beta51*z)+tv52*xx_val(2)*exp(1i*del_beta52*z)+tv53*xx_val(3)*exp(1i*del_beta53*z)+tv54*xx_val(4)*exp(1i*del_beta54*z)+tv56*xx_val(6)*exp(1i*del_beta56*z);...
    tv61*xx_val(1)*exp(1i*del_beta61*z)+tv62*xx_val(2)*exp(1i*del_beta62*z)+tv63*xx_val(3)*exp(1i*del_beta63*z)+tv64*xx_val(4)*exp(1i*del_beta64*z)+tv65*xx_val(5)*exp(1i*del_beta65*z)];

[zz,xa_var]=ode45(f,[0 1],init_condit(nn,:),options); %[LP01 LP11a] represent the initial condition

Unrecognized function or variable 'init_condit'.

[Torsten]

Your code should be execuable and reproduce the problem you are describing. Further - for us uneducated forum people - you will need to explain how energy conservation can be checked from the results you get.

[David Goodmanson]

Hi PS, is it true that (before you start changing things around) all of the del_beta variables are real just like all of the tv variables are, and that for example tv23 = tv32 and del_beta23 = -del_beta32, and the same for the others?

Answer 1

[David Goodmanson]

Hello PS,

Assuming the following (before variables get changed around)

tv_jk = tv_kj,                for all j,k  all real
del_beta_jk = -del_beta_kj    for all j,k  all real

and that for the 6x1 solution vector xx_val = y, what you refer to as power is

power = <y|y> = norm(y)^2  = y'*y = sum(y.*conj(y))    

Then:

The equation is

i*dydt = H*y

and with the assumptions on variables above, H is Hermitian, H = H' and consequently

<y|y> is constant in time.

I am not sure what you mean by 'one of the input parameters' not being real, but if parameters are changed in H to make it not Hermitian, then in general <y|y> will not be conserved. (There may be some odd exception where it's still conserved, possibly). But if you change some of the tv to complex and keep with the condition

tv_jk = conj(tv_kj)     % for all j,k

then H remains Hermitian and <y|y> should be constant in time.