Why am I getting imaginary values with ode45 ??

%This is my program R=1;

L=10^-2;

C=400/62505000010;

Is=10^-12;

Vdo=25*10^-3;

f=@(t,y) [(-R*y(1)/L-(y(2)/(C*L)-(Vdo/L)*log(1+(y(1)/Is))));y(1)];

close all;

tic

options = odeset('RelTol',10^-12,'AbsTol',10^-8);

[t,y]=ode45(f,0:10^-3:0.01,[10^-6 -10^-6],options);

toc

plot(t,y)

xlabel('time "t"')

ylabel('output "y"')

figure(2)

plot(y(:,2),y(:,1),'r')

xlabel('current I');

ylabel('charge Q')

I am simulating a simple circuit, which is governed by two differential equations

   differential (Q)=I;

   differential(I)=(-R/L)*I-(Vdo/L)*log(1+(I/Is))-Q/(C*L);
I,Q are my variables. I am getting the simulated output but the problem is I am getting imaginary values 
for y.
As the only possible cause could be 'log' term in the second differential equation,
 from that         1+I/Is>0     ==> I>-Is   'to get real values'
so I used odeset to put obsolute tolerence limit as 10^-8.
But I am still getting imaginary values for 'y'!!  Why??

Please let me know what I am doing wrong !!

Answer 1

[Shoaibur Rahman]

Your I has both negative and positive values, which are in the range of 10^(-7), but Is is in 10^(-12). So, I/Is gives high negatives for some negative values of I. In that case, log is getting a negative input, and thus resulting imaginary parts.

I is the current, and negative indicates the direction. So, when using inside log, you can use both I and Is in same direction, hence same signs. This is also equivalent to making I as abs(I) inside the log.

Answer 2

[Sai charan Bandi]

Thank You for replying Shoaibur Rahman,

My question how you are saying I is taking negative values also ?? Could you elaborate a little please!! and using abs(I) is the option that I know but I want I to take values > -Is with out using any absolute values just by changing the options in odeset

and

'AbsTol' is not taking negative values, even then when I am using AbsTol as 10^-8 how can I take negative values ??

       as AbsTol > - Is   i.e              10^-8>-10^-12
maybe I am missing some understanding of ode45 here! please explain me if you can understand my question !!

[Shoaibur Rahman]

One column of y is current, right? Look at the plot in figure(2), charge vs. current. You will observe that current has both positive and negative values.

I don't know which circuit you use to develop this system of equations. I guess there may be one or more diode or op-amp elements in the circuit. If that is true, I doubt the log(.) part of your equation would be log|.|, where . indicates the absolute value of .

I am also not an expert user of ode45, but I understand that, in current set up of your system, the only source imaginary parts is log(.)

[Jan]

@Sai charan Bandi: Please use the section for comments and not the one for answers, if you post a comment. Thanks.

[Sai charan Bandi]

Sorry! That's my bad to use answer rather than comment section. @Jan Simon

@Shoaibur Rahman I am using a simple series circuit with R,L,C and one diode without any voltage source.Anyway thanks for the reply.

I could not understand how to control ode45 from not taking any imaginary values.

Answer 3

[Jan]

If you create a function instead of using an anonymous function, you can add a short test for imaginary values:

function dy = fun(t,y)
dy = [(-R*y(1)/L-(y(2)/(C*L)-(Vdo/L)*log(1+(y(1)/Is)))); ...
      y(1)];
if any(~isreal(dy))
  keyboard
end

[Sai charan Bandi]

I know I am getting imaginary values but how to test the performance of ode45 ??

Can we use breakpoints while ode45 executing its integration operation ??

and How can I control ode45 from not taking any imaginary values ??

-Thank You

[Jan]

You cannot tell ODE45 to ibnore imaginary values. When the function to be integrated replies imaginary values, the integration result must also. This is the mathematical definition of an integration.

AbsTol has no relation to the actual values. It means the difference between two different solutions, which are calculated by an order 4 and an order 5 method (see the name ODE 4 5). If the results for one step differs absolutely by less than AbsTol, the step is accepted.

Of course you can use breakpoints - simply try it.

The question about the performance is new, so please ask a new question.

[Sai charan Bandi]

I tried using break points. If I put break point on ode45 line the next step matlab will perform the entire ode calculation and goes to next step. I am asking, is it possible to see , updation of values for each single time-step ??

I cannot tell ode45 to ignore imaginary values , then can I plot the ode output values with imaginary values included . Why my plot function is neglecting the imaginary values of output.??

And finally in my problem I am taking two parameters Q ,I as ode outputs.In One column, I am getting imaginary values and for other column all zeros,which is wrong!! That is why I want to see how ode executing output values in each time step. Is it possible??

-Thanks for the reply

[Shoaibur Rahman]

You can have a look at the ode45 function itself to see how the function is written. Type edit ode45 on your command window, and see that. It is a complex algorithm though.

Regarding the ignoring of imaginary values, the plot function ignores those, not the ode45. This is the way how plot function works. If you are interested to plot the imaginary part then use imag function. However, in RLC circuit analysis, the most effective way will be to see the magnitude and phase, which can be computed using abs and angle functions, respectively.

Two columns in output have some values; none of them is entirely zero, but one is very small comparing to another. This is not surprising since, I = dQ/dt, your first equation. So, the relation between I and Q vastly depends on the time step.