![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1072155/image.png)
solving a differential equation using ode45 but the problem is i didn't get what i expected.
1 view (last 30 days)
Show older comments
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1071950/image.png)
i want to solve this, where σ = 0.1, L0 = 0.5, k= 0.01, tc = 70E-9
and U(Φ0) is
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1071955/image.png)
so that the phase difference Φ0 vs t graph and potential graph U(Φ0) vs Φ0 look like this ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1071960/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1071960/image.png)
I m getting this
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1071965/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1071970/image.png)
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?
0 Comments
Accepted Answer
Sam Chak
on 20 Jul 2022
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1072155/image.png)
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')
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!