solving a differential equation using ode45 but the problem is i didn't get what i expected.

i want to solve this, where σ = 0.1, L0 = 0.5, k= 0.01, tc = 70E-9

and U(Φ0) is

so that the phase difference Φ0 vs t graph and potential graph U(Φ0) vs Φ0 look like this

I m getting this

for Φ0 vs t my program is

=================================================================

ti = 0; % inital time

tf = 10E-5; % final time

tspan = [ti tf];

o =10E5;

k = 0.01; % critical coupling strength

L = 0.5;

s = 0.1;

tc = 70E-9; % photon life time in cavity

f = @(t,y) [(-(s^2)*k/tc)*sin(y - pi/2) + L*(s^2)/(2*tc)*sin(y + pi/2)/sqrt(1 + cos(y + pi/2))];

[T, Y] = ode45(f,tspan, 0);

Y = linspace(-3,3,length(Y));

U = zeros(length(Y),1) ;

for i = 1:length(Y)

U(i) = -(1E-6).*(Y(i)) - (2 + (0.1)^2 ).*((0.033)./(70E-9)).*cos(Y(i) - pi/2) + (0.5).*(((0.1)^2)./(70E-9)).*((1 + cos(Y(i) + pi/2))^(0.5));

end

plot(Y,U);

ax = gca;

ax.XAxisLocation = 'origin';

ax.YAxisLocation = 'origin';

xlabel('\phi')

ylabel('U')

=======================================================================

for U( Φ0) vs Φ0

==========================================================================

ti = 0; % inital time

tf = 10E-5; % final time

tspan = [ti tf];

o =10E5;

k = 0.01; % critical coupling strength

L = 0.5;

s = 0.1;

tc = 70E-9; % photon life time in cavity

f = @(t,y) [(-(s^2)*k/tc)*sin(y - pi/2) + L*(s^2)/(2*tc)*sin(y + pi/2)/sqrt(1 + cos(y + pi/2))];

[T, Y] = ode45(f,tspan, 0);

Y = linspace(-3,3,length(Y));

U = zeros(length(Y),1) ;

for i = 1:length(Y)

U(i) = -(1E6).*(Y(i)) - (2 + (0.1)^2 ).*((0.033)./(70E-9)).*cos(Y(i) - pi/2) + (0.5).*(((0.1)^2)./(70E-9)).*((1 + cos(Y(i) + pi/2))^(0.5));

end

plot(Y,U);

ax = gca;

ax.XAxisLocation = 'origin';

ax.YAxisLocation = 'origin';

xlabel('\phi')

ylabel('U')

====================================================================================

please tell me what's wrong in this program and what parameter can make my result correct?

Answer 1

[Sam Chak]

@SAHIL SAHOO

Nothing wrong with code. The ode45 produces the solution based your given equation and input parameters.

However, one thing is obvious though. This initial value is non-zero.

f = @(t,y) [(-(s^2)*k/tc)*sin(y - pi/2) + L*(s^2)/(2*tc)*sin(y + pi/2)/sqrt(1 + cos(y + pi/2))];
[T, Y] = ode45(f,tspan, 0);
subplot(211)
plot(T, Y), xlabel('t'), ylabel('\sigma')
phase = linspace(-pi, pi, 3601);
U = - (1E-6)*phase - (2 + 0.1^2)*((0.033)/(70E-9))*cos(phase - pi/2) + 0.5*((0.1^2)/(70E-9))*sqrt(1 + cos(phase + pi/2));
subplot(212)
plot(phase, U), xlabel('\sigma'), ylabel('U')

[SAHIL SAHOO]

did i solve the equations correctly ? can you verify that?

[Sam Chak]

Hi @SAHIL SAHOO, Other than your initial value (0) is obviously different, I don't have a reliable source to validate. Maybe the red curves are wrong since they didn't refer to your equation and initial value.

[SAHIL SAHOO]

is there any other method other than ode to solve the equation i mentioned above?

[Sam Chak]

@SAHIL SAHOO, Other methods will still give the same result as ode45 because your initial phase is ZERO radian while the one you compared with is NON-ZERO. Therefore, you can never obtain the result like the red curve.

... unless you ensure that everything (equations and parameters) is absolutely 100% the same.

[SAHIL SAHOO]

okay, thanks for you help man.

[Torsten]

@SAHIL SAHOO

I don't understand why you overwrite your results from ode45 for Y here

[T, Y] = ode45(f,tspan, 0);
Y = linspace(-3,3,length(Y)); 

and here

[T, Y] = ode45(f,tspan, 0);
Y = linspace(-3,3,length(Y)); 

[SAHIL SAHOO]

those are the two different progrsm that i made, one for the Φ0 vs t and second one is for the U( Φ0) vs Φ0

[Torsten]

But you never use the solution for Φ0 from ODE45 in both cases - you overwrite it just after you got it, prescribe an arbitrary range for Φ0 (Y = linspace(-3,3,length(Y))) and calculate U within this range.

[Sam Chak]

Hi @Torsten, I guess @SAHIL SAHOO want to plot U as a function of the phase variable within the range.

[Torsten]

Maybe. But why does he/she solve the ODE ?

[SAHIL SAHOO]

yeah maybe this is the problem, that's why I can't solve the.

can you suggest a method to solve U for the value of Φ0 i obtained from the ODe45?

but what about the Φ0 vs time? is it correct? i change the initial condition still not affecting the outcome.

[SAHIL SAHOO]

@Torsten is right, I want to plot the U as function of Φ0 that i obtained from the ODe45, which is in the range of Φ0 = (-3 ,3)

[Torsten]

can you suggest a method to solve U for the value of Φ0 i obtained from the ODe45?

Read my answer given here:

https://de.mathworks.com/matlabcentral/answers/1764070-is-the-for-loop-is-wrong-what-can-be-the-solution?s_tid=prof_contriblnk

[Sam Chak]

Also this part suggests a piecewise continuous thingy.

But I highly suspect that it has something to do with this term because it happens at .

[SAHIL SAHOO]

if i don't want to use the Φ0 i obtained from the ODe45 to get U then i have to use for loop right?

[Torsten]

No.

Z = linspace(-3,3,size(Y,1)); 
U = -1E6*Z - (2 + 0.1^2 )*0.033/70E-9*cos(Z - pi/2) + 0.5*0.1^2/70E-9*(1 + cos(Z + pi/2)).^0.5;   
plot(Z,U)