I'm trying to plot some level curves from the differential equation dy/dx=-(x^2-x)/(y^2-2y) using dsolve. By hand I get the implicit solution (1/3)*y^3 -y^2 = -(1/3)*x^3+(1/2)*x^2+c, wich not even my HP50g finds an explicit solution.
With MATLAB I find an explicit, complex answer, but that way fplot nor ezplot are able to plot the curves
here's the code I wrote
%%%%
syms y(x) x
eqn = diff(y) == (-x^2+x)/(y^2-2*y)
y0 = [-3 -2 -1 1 2 3]
for k=1:length(y0)
cond = y(0) == y0(k)
sol = dsolve(eqn,cond)
ezplot(sol)
hold on
end
I've been able to plot it with ode15s, but it doesn't give a smooth curve, since it only plots the solution interval containing the initial condition. Also tried plotting fplot(real(sol)), but at the vertical asymptotes, the function looks kind of mirrored.

 Accepted Answer

Try this:
syms y(x) x
eqn = diff(y) == (-x^2+x)/(y^2-2*y);
y0 = [-3 -2 -1 1 2 3];
for k = 1:length(y0)
cond = y(0) == y0(k);
sol{k} = dsolve(eqn,cond);
af{k} = matlabFunction(sol{k});
end
cm = colormap(jet(numel(y0)));
axh = axes('NextPlot','Add');
x = linspace(-2*pi, pi, 150);
for k = 1:numel(af)
fcn = af{k};
plot(x, real(fcn(x)),'-', 'Color',cm(k,:))
plot(x, imag(fcn(x)),'--', 'Color',cm(k,:))
grid
end
It creates anonymous functions from your ‘sol’ results, then uses them to plot the real and imaginary parts in a separate loop.

8 Comments

Thanks man, this really improves the model. But that doesn't quite get the hole function. Given that there are 3 solution intervals, and it does a great job on the x axis, it still doesn't account for the other vertical asymptote (or others on higher degree functions). Also, without previously knowing the shape of the level curves, it's still very confusing, right? What i really think could solve it was to get an implicit solution from dsolve, but after a couple os days of research, didn't find anything like that, wich makes me wonder if such a thing even exists... Do you know something about that? Here's what it should look like, plotting the implicit function on DESMOS
My pleasure.
I have no idea what you are doing. The ‘image.jpg’ you posted does not appear to differentiate the real and imaginary solutions.
My goal is to get your code to produce and plot the expected real and imaginary results from it, and I have.
If your equations are functions of two variables and produce a result that has a third dimension (essentially a parametric ‘z=f(x,y)’ where both ‘x’ and ‘y’ are parametric functions) code it that way and use the contour function to plot the isolevels.
Yes, your plot does exactly what you said. I'm sorry, the problem might not be perfectly clear. My goal is to plot the solution function from the ODE with different initial conditions, wich should look like as the image I've uploded. With ode15s i've been able to get the general idea about the behavior from the ODE solutions, but i can't get the continuous functions that are expected (yes, they're only continuous by parts, but with the implicit function it gets the smooth curves I need)
Here's the ode15s code
if true
% code
syms y(x) x
eqn = @(x,y)(-x^2+x)/(y^2-2*y)
y0 = -2:0.25:5
for k=1:length(y0)
a=5;
for z=0:1
tspan1 = [z a];
tspan2 = [z -a];
[X1,Y1] = ode15s(eqn, tspan1, y0(k));
[X2,Y2] = ode15s(eqn, tspan2, y0(k));
plot(X1,Y1)
hold on
plot(X2,Y2)
axis([-3 3 -1 3])
end
end
end
As the function isn't symetrical with respect to the x axis, I don't have a continuous curve at the vertical asymptotes. That's when I tried the 'dsolve' solver, and it gave me the explicit solutions, with the complex functions that you already solved. What I wanted though, as an output from 'dsolve' was either a continuous explicit function (wich is obviously not the case), or the implicit function, so I can plot it with various initial conditions. Thank you again.
I am lost.
That seems to be the result you want.
What else needs to be done?
I have the ODE, defined earlier. It's a simple separable equation, which solved by hand gives the implicit solution (1/3)*y^3 -y^2 = -(1/3)*x^3+(1/2)*x^2+c. Let's say I need to plot the solution with initial conditions y(0)=1. Solving for c gives c=-2/3, and I get the following plot:
This is the whole function, though the EDO solution lies only in the interval of y from 0 to 2.
with ode15s I get the following:
which is the correct solution, but limitates me on the interval of solution. To get the whole function, I should find the other initial conditions, and that was never my idea.
if true
% code
syms y(x) x
eqn = @(x,y)(-x^2+x)/(y^2-2*y)
y0 = [-0.729 1 2.735];
a=5;
z=0;
for k=1:length(y0)
tspan1 = [z a];
tspan2 = [z -a];
[X1,Y1] = ode15s(eqn, tspan1, y0(k));
[X2,Y2] = ode15s(eqn, tspan2, y0(k));
plot(X1,Y1)
hold on
plot(X2,Y2)
axis([-2 2.5 -2.5 4])
end
end
As you can see, by forcing initial conditios that I "shouldn't know", not without the implicit plot. Finally, with dsolve, I get
if true
% code
syms y(x) x
eqn = diff(y) == (-x^2+x)/(y^2-2*y)
y0 = [1]
for k=1:length(y0)
cond = y(0) == y0(k)
sol = dsolve(eqn,cond)
ezplot(sol)
axis([-2 2.5 -2.5 4])
hold on
end
end
which I know is correct, but again, only gives me the solution on the defined interval of y.
What I was wondering is, if there's any way to get from dsolve an implicit function whenever there's no explicit continuous solution, so I could plot the curves such as the first image on this comment.
I believe dsolve has done everything it is capable of doing. Other symbolic engines, such as Maple, may be able to provide more solutions.
Ok, I can live with the numerical solution, I guess. Thanks!
As always, my pleasure.
If my Answer helped you solve your problem, please Accept it!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!