Getting complex solution with Matlab ode45
Mostrar comentarios más antiguos
I'm trying to solve such an ODE:
with boundary condition
I have such a problem that if the solution of y would come with complex part such as 1.1462 + 0.0000i. However, if I add an abs in sqrt there will be no longer complex parts, but the solution is not what I desires. Is there anything I can do to avoid this?
My matlab code is
clearvars;
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
tspan2=[a 0.99];
options=odeset('Reltol',1e-10);
[t,y]=ode45(odefun_mon_pc1,tspan2,0.0001,options);
h_pc_mon=plot(t,y,'--blue','linewidth',1);
10 comentarios
Paul
el 6 de En. de 2024
What is the desired solution?
Kaiwei
el 6 de En. de 2024
Kaiwei
el 6 de En. de 2024
Let's check
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
%tspan2=[a 0.99];
%options=odeset('Reltol',1e-10);
%[t,y]=ode45(odefun_mon_pc1,tspan2,0.0001,options);
%h_pc_mon=plot(t,y,'--blue','linewidth',1);
The derivative at t = a and y = 0 is
odefun_mon_pc1(a,0)
But, y just a bit more than 0, like the initial condition in the code
odefun_mon_pc1(a,0.0001)
So it doesn't take much to make the derivative negative.
Kaiwei
el 6 de En. de 2024
Kaiwei
el 6 de En. de 2024
Kaiwei
el 6 de En. de 2024
Respuestas (2)
It looks like your differential equation is very sensitive (I mean that loosely)
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
tspan2=[a 0.99];
Here's the initial code:
options=odeset('Reltol',1e-10);
[t,y] = ode45(odefun_mon_pc1,tspan2,0.0001,options);
figure
h_pc_mon=plot(t,y,'--blue','linewidth',1);
Use the stated boundary condition:
options=odeset('Reltol',1e-10);
[t,y] = ode45(odefun_mon_pc1,tspan2,0*0.0001,options);
figure
h_pc_mon=plot(t,y,'--blue','linewidth',1);
Use the stated boundary condition with a smaller step size
options=odeset('Reltol',1e-10,'MaxStep',1e-6);
[t,y] = ode45(odefun_mon_pc1,tspan2,0*0.0001,options);
figure
h_pc_mon=plot(t,y,'--blue','linewidth',1);
If you cut tspan down to focus on the what's happening for t < 0.3, you'll probably find some interesting behavior.
I didn't try setting the MaxInitialStep to a small value; that's usually a good idea. And then maybe the MaxStep won't have to be set so small. You'll have to do some more analysis.
15 comentarios
Kaiwei
el 6 de En. de 2024
Walter Roberson
el 7 de En. de 2024
Editada: Walter Roberson
el 7 de En. de 2024
format long g
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
tspan2=[a 0.99];
options=odeset('Reltol',1e-10);
[t,y]=ode45(odefun_mon_pc1,tspan2,0.0001,options);
subplot(2,1,1)
plot(t,real(y),'--blue'); title('real');
subplot(2,1,2)
plot(t,imag(y),'--red'); title('imag')
[min(imag(y)), max(imag(y))]
With the default "format short", values in that range show up as 0.0000i
Maybe of interest: For given value of a, y0 must be chosen smaller or equal than the curve value in order to get an increasing solution for the differential equation.
syms y t real
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
sol = solve(t*D*y+A*t+C+sqrt(B*y)==0,'ReturnConditions',1)
sol.y
sol.conditions
y = sol.y(1)
tnum = a:0.01:0.99;
format long
ynum = arrayfun(@(tnum)double(subs(y,t,tnum)),tnum);
ynum= ynum(:)
plot(tnum,ynum)
xlabel('tstart')
ylabel('Frontier for y0')
For solutions with absolute value <1, it is important to restrict AbsTol, not RelTol in the computations:
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
tspan2=[a 0.99];
options=odeset('AbsTol',1e-10);
[t,y]=ode45(odefun_mon_pc1,tspan2,0,options);
subplot(2,1,1)
plot(t,real(y),'--blue'); title('real');
subplot(2,1,2)
plot(t,imag(y),'--red'); title('imag')
[min(imag(y)), max(imag(y))]
In response to this comment, I don't see the solution have any element with a non-zero imaginary part.
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=0.2778;
odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
tspan2=[a 0.99];
options=odeset('Reltol',1e-10,'MaxStep',1e-6);
[t,y] = ode45(odefun_mon_pc1,tspan2,0,options);
any(imag(y))
In response to this comment
sol = ode45(odefun_mon_pc1,tspan2,0,options);
f = @(t) deval(sol,t);
tval = 0.3:.01:.4;
[yval,yprime] = f(tval);
figure
plot(tval,yval,tval,yprime)
legend('y','yprime','Location','Northwest')
Kaiwei
el 7 de En. de 2024
Kaiwei
el 7 de En. de 2024
Actually even setting Abstol you get the same result, the y still contains xx+0.000i part. I'm a little concerned that the appearance of 0.000i is a bug, there must be some cases there is negative numbers in sqrt.
I think once your function returned complex values (maybe in some intermediate step of the integration), the complex part in the solution remains - even if the imaginary part of the integration step didn't contribute to the final solution.
Do you know how to solve an ODE which can not be explicitly expressed as y'=@(t,y) ? (i.e. I want to squre on both side to get rid of the sqrt part. But after doing so I will get an ODE expressed as f(t,y,y')=0 ) I don't know how to solve that. Maybe you know how to do that?
Solve
z1' = z2
0 = f(t,z1,z2)
using ode15s with z1 = y and z2 = y'. Use the mass matrix option to do so.
But note that this squaring will complicate things for the solver because you will have to tell him which path to take when solving x^2 = c and c becomes 0 - the path of the positive or the negative square root.
Paul
el 7 de En. de 2024
Please provide a specific example that illustrates the "appearance of 0.000i". As shown, when the code runs here on Answers as in the first part of this comment the imaginary part of all elements of y is identically zero.
Kaiwei
el 7 de En. de 2024
Paul
el 7 de En. de 2024
Without seeing an actual example what you're seeing that "appears," I'm afraid I can't be of any further help. Maybe someone else can. Why not post the actual code and output so we all know what you're talking about?
Kaiwei
el 7 de En. de 2024
Here is the "protocol" of the computation:
format long
main()
function main
A=-4.32;
B=4;
C=1.2;
D=0.8;
a=5/18;
%odefun_mon_pc1=@(t,y) (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D);
tspan2=[a 0.99];
options=odeset('AbsTol',1e-10);
[t,y]=ode45(@odefun_mon_pc1,tspan2,0,options);
function dy = odefun_mon_pc1(t,y)
y
dy = (t*D*y+A*t+C+sqrt(B*y))/((t*t-t)*D)
end
end
Yes, I thought we already established that original code eventually leads to an integration step where y becomes negative and dy becomes complex valued, at which point y will have a non-zero imaginary component. So I'm still not sure what the concern is with the original code.
Maybe you're just referring to the display of the output in the Matlab window, which depends on the format of the output display and how Matlab diplays 0 values (for either real or imaginary parts) of complex-valued arrays.
Consider the real array
x = [0;eps]
Notice that with the defaut formating the first element is displayed as "0", which indicates that it's exactly zero ( and multiplied by 1e-15)
But, with a complex array
y = [0; 0 + 1i*eps; eps]
Ther real and imaginary parts of elements where those parts are exactly zero are displayed with "0.0000".
But if we change format, we see those parts are, in fact, exactly zero
format long e
y
Hi @Kaiwei
As noted by other users, the complex-valued solution is caused by the square root term 'sqrt(B*y)' when the output y is negative. In theory, the integration will yield a negative output y when y-prime (
) is negative. However, the output y vector in the numerical solution does not have a negative real part. Thus, in a human sense, the argument in the square root term should remain positive, as you would expect a monotonically increasing pattern in the output.
) is negative. However, the output y vector in the numerical solution does not have a negative real part. Thus, in a human sense, the argument in the square root term should remain positive, as you would expect a monotonically increasing pattern in the output.Here, using the conventional Runge-Kutta 4th-order (RK4) solver, we can examine what really happens to the initial step when the step size is very small. The solution of RK4 consists of 4 terms, 
, 
, 
, and 
. The result below shows that 
becomes complex-valued because become negative.
, , , and . The result below shows that becomes complex-valued because become negative.format long g
%% Kaiwei's model
% Parameters
A = -4.32;
B = 4;
C = 1.2;
D = 0.8;
a = 5/18; % 0.277777777777777777777777777777777
% Ordinary Differential Equations
odefcn = @(t, y) (t*D*y + A*t + C + sqrt(B*y))/(D*(t.^2 - t));
% Simulation time span
tStart = a; % 5/18
h = 0.00000013/18; % a very small step size
tEnd = a + h; % check what happen to the initial step
% tEnd = 1 - 0.0013/18; % singularity occurs at 18/18 = 1
t = tStart:h:tEnd;
% Initial value
y0 = 0;
% Call RK4 Solver
yRK4 = RK4Solver(odefcn, t, y0);
% Plotting the solutions
plot(t, yRK4, '-o')
grid on
xlabel('\alpha')
ylabel('U')
%% Runge-Kutta 4th-order Solver
function y = RK4Solver(f, x, y0)
y(:, 1) = y0; % initial condition
h = x(2) - x(1); % step size
n = length(x); % number of steps
for j = 1 : n-1
k1 = f(x(j), y(:, j))
k2 = f(x(j) + h/2, y(:, j) + h/2*k1)
k3 = f(x(j) + h/2, y(:, j) + h/2*k2)
k4 = f(x(j) + h, y(:, j) + h*k3)
y(:, j+1) = y(:, j) + h/6.0*(k1 + 2*k2 + 2*k3 + k4)
end
end
3 comentarios
Kaiwei
el 7 de En. de 2024
Hi @Kaiwei
An implicit ODE doesn't look like yours. Essentially, as @Torsten explained, once the complex-valued solution appears in y, the imaginary part 'cannot be removed'. The values of the imaginary part do not contribute significantly because they are very, very small. You can evaluate the values of the imaginary part below:
By the way, just out of curiosity, what physical phenomenon does your differential equation describe?
format long g
%% Kaiwei's model
% Parameters
A = -4.32;
B = 4;
C = 1.2;
D = 0.8;
a = 5/18; % 0.277777777777777777777777777777777
% Ordinary Differential Equations
odefun = @(t, y) (t*D*y + A*t + C + sqrt(B*y))/(D*(t.^2 - t));
% Simulation time span
tStart = a; % 5/18
h = 0.13/18; % a very small step size
tEnd = 1 - 0.13/18; % singularity occurs at 18/18 = 1
t = tStart:h:tEnd;
% Call ode45 solver
tspan = a:0.13/18:(1 - 0.13/18);
y0 = 0;
options = odeset('Reltol', 1e-10, 'MaxStep', (13/18)/1e3);
sol = ode45(odefun, tspan, y0, options);
[y, yp] = deval(sol, tspan) % yp is y-prime (y')
plot(tspan, y, tspan, yp), grid
xlabel('\alpha'), legend('y', 'y-prime')
Kaiwei
el 7 de En. de 2024
Categorías
Más información sobre Numeric Solvers en Centro de ayuda y File Exchange.
Productos
el 6 de En. de 2024
el 7 de En. de 2024
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
