Plot differential equations with respect to two variables in a 3d plane

I have a system of three differential equations and coded as follows.

I have a seperate function "first_term.m" to create the first part of the equation and another function "second_term.m" to create the second part. And then there is another function "add_RHS.m" to combine both these terms and pass it to ModelRHS(t,x,param).

Here's my add_RHS.m function that defines model equations.

if some condition > 0
    dxdt(j) = dxdt(j) + first_term ;
end

if some condition < 0
    dxdt(j) = dxdt(j) + second_term ;
end

Both these first and second terms consist of and . I need to plot vs vs in a 3D plane. Can someone plese suggest a way to do this? Simulation for this system of differential equations is given below.

editparams; %file that contains parameters
Tend = 100;
Nt = 100;
% Define RHS functions
RHS = @(t,x)RHS(t,x,param); 
%Execution-----------------------------------------------------------------
x0 = [0.004, 0.05, 0.1]; %Initial condition
t = linspace(0,Tend,Nt); %TSPAN
[t, A] = ode45(RHS, t, x0);

[Sam Chak]

@UserCJ, Are you looking for this?

plot3(X, Y, Z)

https://www.mathworks.com/help/matlab/ref/plot3.html

[UserCJ]

No, I'm looking more likely for a surface plot.

Answer 1

[Davide Masiello]

You can use surf

https://ch.mathworks.com/help/matlab/ref/surf.html

or contourf

https://ch.mathworks.com/help/matlab/ref/contourf.html

[UserCJ]

@Davide Masiello Thanks for your answer. That means if I need to plot against and I need to consider

A(:,3) 

and plot against other two variables and using surf or contourf ?

My guess is A gives all the x values and I'm having a difficulty of extracting .

[Davide Masiello]

You can evaluate the derivative using deval

https://ch.mathworks.com/help/matlab/ref/deval.html

[UserCJ]

Thanks @Davide Masiello here's what I did

A= [t, sol];
tvec = linspace(0,Tend, 100)
[sxint,spxint]=deval(A,tvec);

and I also tried

[t, sol] = ode45(RHS, t, x0);
[sxint,spxint]=deval(sol,tvec);

but in both cases I got the same error message as

Error using deval (line 46)
sol must be a structure returned by a differential equation solver.

Can you please tell me what I am missing?

[Torsten]

sol = ode45(RHS, t, x0);
tvec = linspace(0,Tend, 100)
[y,yp] = deval(sol,tvec);
x1 = y(:,1);
x2 = y(:,2);
dy3 = yp(:,3);
plot3(x1,x2,dy3)

Using "surf" or "contour" does not make sense here.

[UserCJ]

@Torsten Thanks, so much for your answer. I'd like to get this plot as a surface plot. How can I make it?

[Torsten]

What kind of surface do you expect to get if you only have a curve (x1(t),x2(t)) in the x1-x2-plane ?

Is it this what you mean ?

https://de.mathworks.com/matlabcentral/answers/368601-plot-a-curve-as-a-vertical-3d-surface-in-matlab

[UserCJ]

@Torsten yes, that's what I'm looking for!

[Torsten]

Then use Star's code supplied there.

[UserCJ]

@Torsten but still it has only one euation and not like mine.

[Torsten]

x = linspace(0, 2*pi, 50);                              % Independent Variable
y = sin(x);                                             % Dependent Variable
z = [(1:-0.2:0)'*ones(size(x))];                        % ‘Z’ Matrix

in Star's code.

x is x1 in your case, y is x2 in your case, an appropriately modified z is dy3 in your case.

Can you take it from here ?

[UserCJ]

@Torsten Thanks so much for your help! Does this makeany sense ? I'm so sorry I'm new to this and still learning.

I realised that , , takes only three values in each. From what I've understood, I need to make them as matrices?

tvec = linspace(0,Tend, 100);
[y,yp] = deval(sol,tvec);
x1 = y(:,1);
x2 = y(:,2);
dy3 = yp(:,3);
surf(repmat(x1,size(dy3,1),1), repmat(x2,size(dy3,1),1), dy3)
grid on
axis equal

[Torsten]

sol = ode45(RHS, t, x0);
tvec = linspace(0,Tend, 100)
[y,yp] = deval(sol,tvec);
x1 = y(:,1).';
x2 = y(:,2).';
dy3 = yp(:,3).';
nsteps = 20;
DY3 = zeros(nsteps,numel(x1));
for i=1:numel(x1)
    DY3(:,i) = linspace(0,dy3(i),nsteps).'
end
figure(1)
surf(repmat(x1,size(DY3,1),1),repmat(x2,size(DY3,1),1),DY3)
grid on
axis equal

[UserCJ]

@Torsten Thanks so much for your help. I tried according to your code and heres what I got.

However if I just run the first run the code until for loop and then try surf in the following way I got a different image. But I think this time the axes are different.

surf(DY3)
shading interp

Btw,

surf(DY3)

does not work together with the first half of the code (until the end of for loop)

Do you haveany idea of what's happening here?

[Torsten]

I don't have

surf(DY3)

in the code I provided.

Further, I requested 100 data points from ODE45 by the commands

tvec = linspace(0,Tend, 100)
[y,yp] = deval(sol,tvec);

In your plot, I see only three.

The second plot cannot be what you want. The surface plot should be a curve in the x1x2 plane pulled in the dy3 direction.

[UserCJ]

@Torsten Here's the code I used! I changed only the "nsteps = 100;" from yours. And I'm not sure where that 3 is coming from.

tvec = linspace(0,Tend, 100);
[y,yp] = deval(sol,tvec);
x1 = y(:,1).';
x2 = y(:,2).';
dy3 = yp(:,3).';
nsteps = 100;
DY3 = zeros(nsteps,numel(x1));
for i=1:numel(x1)
    DY3(:,i) = linspace(0,dy3(i),nsteps).';
end
figure(1)
surf(repmat(x1,size(DY3,1),1),repmat(x2,size(DY3,1),1),DY3)
shading interp
grid on
axis equal

[Torsten]

In your first plot, it looks as if the surface is only made up of three data points: one in the neighbourhood of (0,0), one in the neighbourhood of (0.2,0.2) and the third nearby.

And in a previuos post you reported that x1,x2 and dy3 have only three values each.

So this is not the case here when you try to plot the surface ? x1,x2 and dy3 are row vector with 100 elements each ?

[UserCJ]

@Torsten Let me get this thing from the beginning;

size([y,yp]) = 3     200 %since Tend= 200
size(x1) = 1    3
size(x2) = 1    3
size(dy3) = 1     3
size(DY3) = 100   3

Where am I going wrong with three points? How can I get x1,x2 and dy3 as row vectors with 100 elements each?

[Torsten]

Use

tvec = linspace(0,Tend,100).'
[y,yp] = deval(tvec,sol);

y and yp should be 100x3.

Is this correct ?

[UserCJ]

@Torsten yes, that's correct

size(y) = 3   100
size(yp) = 3   100

Still I get the same output.

[Torsten]

No, it's not correct - it should be 100x3.

But if deval works like this, change the code to

tvec = linspace(0,Tend, 100);
[y,yp] = deval(sol,tvec);
x1 = y(1,:);
x2 = y(2,:);
dy3 = yp(3,:);
nsteps = 10;
DY3 = zeros(nsteps,numel(x1));
for i=1:numel(x1)
    DY3(:,i) = linspace(0,dy3(i),nsteps).';
end
figure(1)
surf(repmat(x1,size(DY3,1),1),repmat(x2,size(DY3,1),1),DY3)
shading interp
grid on
axis equal

[UserCJ]

@Torsten I changed the code and here's the output

Still

size(y) = 3  100
size(yp)= 3    100

[UserCJ]

@Torsten Thank you so much for your help!

[Torsten]

Looks better, but maybe you could scale dy3 by a factor of 10000.

[UserCJ]

Thanks @Torsten. Yes, I did. It took so long to get in here!

Answer 2

[Bruno Luong]

Try this (adapt to your code):

%editparams; %file that contains parameters
Tend = 10;
Nt = 100;
% Define RHS functions
RHS = @(t,x) sin(t).*x.^2;
%Execution-----------------------------------------------------------------
x0 = rand; %Initial condition
t = linspace(0,Tend,Nt); %TSPAN
[t, x] = ode45(RHS, t, x0);

close all
tgrid = t;
Nx = 60;
xgrid = linspace(min(x),max(x),Nx);
[T,X] = meshgrid(tgrid,xgrid);
DXDT = RHS(T,X);
surf(tgrid,xgrid,DXDT);
hold on
dxdt = RHS(t,x);
plot3(t,x,dxdt,'r','Linewidth',3);

[UserCJ]

@Bruno Luong Thanks, so much for this. It returned the following error;

Error using surf (line 71)
Z must be a matrix, not a scalar or vector.

[Bruno Luong]

My latest code runs on TMW onlne server (see the plot)

[UserCJ]

@Bruno Luong Thanks for your help! This is how I did in my code.

[t, sol] = ode45(RHS, t, x0, options);
%A= [t, sol];
tgrid = t;
Nx = 100;
xgrid = linspace(min(sol(:)),max(sol(:)),Nx);
[T,X] = meshgrid(tgrid,xgrid);
DXDT = RHS(T,X);
surf(tgrid,xgrid,DXDT)

When I checked every single line, I can data matrix for xgrid, but then DXDT returns only one set of values.

0.0538
0.0168
0.0105

I think that is why it's returning an error. Do you have any idea how I could fix this?Thanks in advance!

[Bruno Luong]

Change your function RHS so it can deal with multiple inputs.

You haven't provided your RHS code so nobody can help you on the dark side.

Another way is to change the statement DXDT = ... by:

DXDT = arrayfun(RHS,T,X);

[Torsten]

It won't work in your case since you solve for 3 unknown functions, not 1 as in Bruno's example code.

[Bruno Luong]

In reply to your code

editparams
if some type of edge
if another type of edge        
introduce parameters
addRHS

I would said modify my code to

Do something with your RHS to compute DXDT correctly