ODE45 solver problem outcome

I have 5 differential equations and a couple of extra equations that I want to solve for 5 variables.

The problem is that if you look at the graph you see that the outcomes do not match. dh/dt should be -u1(t) but that is not the case in the graphs. I do not know where it goes wrong in my code.

If I put the equations for u1(t), u2(t) and ud(t) in the Eqns the same results are obtained which seems strange to me. How do I write the code such that all equations are true?

clear all;

ds = 18.32
ws = 45
br = 0.0015
bd = 0.0015
rho = 1019
g = 9.81
us = 0.06
massvessel = ds * ws * rho
wd = 1.13

syms u1(t) u2(t) ud(t) h(t) wr(t) T Y

u1(t) == (-us * ds + wr(t) * u2(t))/wr(t);
u2(t) == (us*ds + wr(t)* u1(t))/wr(t);
ud(t) == (u2(t)*wr(t))/wd;

Eqns = [diff(u1(t),t) == g*(h(t)/ds) - br * (u1(t)); 
        diff(u2(t),t) == -g * (h(t)/ds) - br * (u2(t));
        diff(ud(t),t) == -g * (h(t)/ws) - bd * (ud(t));
        diff(h(t),t) == - u1(t);
        diff(wr(t),t) == -us];

[DEsys,Subs] = odeToVectorField(Eqns);

DEfcn = matlabFunction(DEsys, 'Vars',{T,Y});
tspan = linspace(0, 200, 251);
Y0 = [0; 0; 0; 0; 30]+0.001;
[T, Y] = ode45(DEfcn, tspan, Y0);

[darova]

Here is the mistake

You are not assigning equation to variable. Use '=' sign once

[Iris Heemskerk]

@darova Thank you so much for your help!

If I change this however, the following error appears:

Error using mupadengine/feval_internal (line 172)
Unable to convert the initial value problem to an equivalent dynamical system.
Either the differential equations cannot be solved for the highest derivatives or
inappropriate initial conditions were specified.

Error in odeToVectorField>mupadOdeToVectorField (line 171)
T = feval_internal(symengine,'symobj::odeToVectorField',sys,x,stringInput);

Error in odeToVectorField (line 119)
sol = mupadOdeToVectorField(varargin);

Error in simulatieexp1nonlinear (line 33)
[DEsys,Subs] = odeToVectorField(Eqns);

Answer 1

[Star Strider]

The ‘==’ are correct here. They are symbolic equations, not logical operations. Otherwise, I would agree.

Also, your code works correctly. The integrated value ‘h(t)’ will appear different from ‘-u1(t)’ because it is integrated. If you plot ‘u1(t)’ (integrated as ‘Y(:,2)’) against the negative of the derivative of ‘h(t)’:

figure
plot(T,Y(:,2), '-b',  T,-gradient(Y(:,4),tspan(2)), '--r')
grid

(with the derivative calculated by the gradient function), they are almost exactly the same.

.

[Iris Heemskerk]

That is weird because the orange line in the graph is the negative gradient of u1(t).

Besides, isn't Y(:,2) u2(t) and not u1(t)? Shouldn't this be:

figure
plot(T,Y(:,1), '-b',  T,-gradient(Y(:,4),tspan(2)), '--r')
grid

Another thing that bothers me is when I change the sequence of the differential equations, the graphs are not the same.

Do you know why this is?

[Star Strider]

‘That is weird because the orange line in the graph is the negative gradient of u1(t).’

That is the purpose of that demonstration. The code does exactly what it should.

‘Besides, isn't Y(:,2) u2(t) and not u1(t)?’

No. That is where the ‘Subs’ output helps:

Subs =
 u2
 u1
 ud
  h
 wr 

That corresponds to the order of the ‘Y’ vector.

Plotting all of them and labeling each one:

figure
for k = 1:size(Y,2)
    subplot(size(Y,2),1,k)
    plot(T, Y(:,k))
    grid
    title(string(Subs(k)))
end

demonstrates this.

‘Another thing that bothers me is when I change the sequence of the differential equations, the graphs are not the same.

Do you know why this is?’

I am not certain what you are doing. However, it is likely that the OdeToVectorField function ls rearranging the variable assignments. Again, look at the ‘Subs’ output (or run my plotting loop). That will tell you how ‘Y’ corresponds to the variables.

.

[Iris Heemskerk]

If i change the code such that only the sequence of the differential equations differs, the outcomes of the graph are totally different. I used the code you wrote for plotting the graphs.

syms wr(t) h(t) u1(t) u2(t) ud(t) T Y

u1(t) == (-us * ds + wr(t) * u2(t))/wr(t);
u2(t) == (us*ds + wr(t)* u1(t))/wr(t);
ud(t) == -(u2(t)*wr(t))/wd;




Eqns =  [diff(wr(t),t) ==  -us;
        diff(h(t),t) == - u1(t);
        diff(u1(t),t) == g*(h(t)/ds) - br * (u1(t)); 
        diff(u2(t),t) == -g * (h(t)/ds) - br * (u2(t));
        diff(ud(t),t) == -g * (h(t)/ws) - bd * (ud(t))]

    DEsys,Subs] = odeToVectorField(Eqns);

DEfcn = matlabFunction(DEsys, 'Vars',{T,Y});

tspan = linspace(0, 200, 251);

Y0 = [30; 0; 0; 0; 0]+0.001;

[T, Y] = ode45(DEfcn, tspan, Y0);

%or this sequence:

syms u1(t) u2(t) ud(t) wr(t) h(t) T Y

u1(t) == (-us * ds + wr(t) * u2(t))/wr(t);
u2(t) == (us*ds + wr(t)* u1(t))/wr(t);
ud(t) == -(u2(t)*wr(t))/wd;




Eqns =  [diff(u1(t),t) == g*(h(t)/ds) - br * (u1(t)); 
        diff(u2(t),t) == -g * (h(t)/ds) - br * (u2(t));
        diff(ud(t),t) == -g * (h(t)/ws) - bd * (ud(t));
        diff(wr(t),t) ==  -us;
        diff(h(t),t) == - u1(t)];

    DEsys,Subs] = odeToVectorField(Eqns);

DEfcn = matlabFunction(DEsys, 'Vars',{T,Y});

tspan = linspace(0, 200, 251);

Y0 = [0; 0; 0; 30; 0]+0.001;

[T, Y] = ode45(DEfcn, tspan, Y0);

I think it is really strange that the order of magnitde on the y-axis is so different in the two graphs.

[Star Strider]

Do the same experiment, however using my loop code:

figure
for k = 1:size(Y,2)
    subplot(size(Y,2),1,k)
    plot(T, Y(:,k))
    grid
    title(string(Subs(k)))
end

to plot them. Note that it uses the ‘Subs’ output to title each plot appropriately.

So the plots for the first sequence:

and the plots for the second sequence:

Only the order is different. Everything else is unchanged.

.

[Iris Heemskerk]

The values on the y-axis are 10000 times larger in the first graph compared to the second one? Furthermore, he initial value for wr(t) does not suffice in the first graph?

[Star Strider]

All I did was to copy and past the code you posted in your previous Comment, then ran it with my plotting loop to create those figures. You apparently changed ‘tspan’ between the code that produced your figures and the code you most recently posted.

The ‘tspan’ that you posted and that I used is:

tspan = linspace(0, 200, 251);

I made no changes in the code you posted. All I did was run both of them as posted, using my subplot loop to plot them (since it titles the subplot plots using the ‘Subs’ output). I did not look at how the ‘Y0’ vector elements corresponded to the order of the function arguments. I left that you you.

.

[Iris Heemskerk]

This is exacty my point. The two codes I posted are the same, the only thing that differs between the two codes is the order of the differential equations. This change in order produces different graphs, while the rest of the code is the same. If you take a close look at the two figures you produced using my code from the previous comment, you see that the graphs are different while the code is the same (except for the order of differential equations).

My only question is, how can the order of differential equations influence the outcome of the solver.

[Star Strider]

That has to do with the way ode45 integrates them. It uses the existing values of the variables in each equation. That value is going to be the result of the previous integration iteration (or the initial conditions in the first iteration), or the newly-evaluated result in the current integration iteration. So the order is significant.

[Iris Heemskerk]

Thank you so much for your help! I suppose that the correct order to use depends on the dependency of the variables of each other and that you should consider which one should be where based on the mathematics.

[Star Strider]

As always, my pleasure!

Exactly!

[Iris Heemskerk]

When one problem is solved another one emerges..

I noticed that there is no difference in graphs if you do not take into account the extra equations. This means that the solver does not use the three additional equations (u1(t), u2(t), ud(t))? I checked this by plotting the ud(t) from the solver against the ud(t) = u2(t)*wr(t))/wd, and the two graphs are not the same. Furthermore, if you put a % before the three additional equations, the outcomes of the graphs are exactly the same.

close all;
clear all;

ds = 18.32
ws = 45
br = 0.0015
bd = 0.0015
rho = 1019
g = 9.81
us = 0.06
massvessel = ds * ws * rho
wd = 1.13

syms h(t) u1(t) u2(t) wr(t) ud(t) T Y

% these 3 equations do not suffice in the graph:
u1(t) == (-us * ds + wr(t) * u2(t))/wr(t);
u2(t) == (us*ds + wr(t)* u1(t))/wr(t);
ud(t) == (u2(t)*wr(t))/wd;

Eqns = [diff(h(t),t) == - u1(t)
        diff(u1(t),t) == g*(h(t)/ds) - br * (u1(t)); 
        diff(u2(t),t) == -g * (h(t)/ds) - br * (u2(t));
        diff(ud(t),t) == -g * (h(t)/ws) - bd * (ud(t));
        diff(wr(t),t) == -us];
        
[DEsys,Subs] = odeToVectorField(Eqns);

DEfcn = matlabFunction(DEsys, 'Vars',{T,Y});

tspan = linspace(0, 200, 251);

Y0 = [0; 0; 0; 0; 30]+0.001;

[T, Y] = ode45(DEfcn, tspan, Y0);

figure
for k = 1:size(Y,2)
    subplot(size(Y,2),1,k)
    plot(T, Y(:,k))
    grid
    legend(string(Subs(k)))
end

figure
plot(T, Y(:,4))
hold on
plot(T, Y(:,3).*Y(:,5)./wd)
legend('ud', 'ud from u2*wr/wd')

[Star Strider]

I dio not understand the problem.

[Iris Heemskerk]

The two codes below produce exactly the same outcome.

This means that the 3 equations are not taken into account by the code.

close all;
clear all;
ds = 18.32
ws = 45
br = 0.0015
bd = 0.0015
rho = 1019
g = 9.81
us = 0.06
massvessel = ds * ws * rho
wd = 1.13
syms h(t) u1(t) u2(t) wr(t) ud(t) T Y

u1(t) == (-us * ds + wr(t) * u2(t))/wr(t);
u2(t) == (us*ds + wr(t)* u1(t))/wr(t);
ud(t) == (u2(t)*wr(t))/wd;

Eqns = [diff(h(t),t) == - u1(t)
        diff(u1(t),t) == g*(h(t)/ds) - br * (u1(t)); 
        diff(u2(t),t) == -g * (h(t)/ds) - br * (u2(t));
        diff(ud(t),t) == -g * (h(t)/ws) - bd * (ud(t));
        diff(wr(t),t) == -us];
        
[DEsys,Subs] = odeToVectorField(Eqns);
DEfcn = matlabFunction(DEsys, 'Vars',{T,Y});
tspan = linspace(0, 200, 251);
Y0 = [0; 0; 0; 0; 30]+0.001;
[T, Y] = ode45(DEfcn, tspan, Y0);

figure

for k = 1:size(Y,2)
    subplot(size(Y,2),1,k)
    plot(T, Y(:,k))
    grid
    legend(string(Subs(k)))
end



close all;
clear all;

ds = 18.32
ws = 45
br = 0.0015
bd = 0.0015
rho = 1019
g = 9.81
us = 0.06
massvessel = ds * ws * rho
wd = 1.13
syms h(t) u1(t) u2(t) wr(t) ud(t) T Y


% these 3 equations do not suffice in the graph:
% u1(t) == (-us * ds + wr(t) * u2(t))/wr(t);
% u2(t) == (us*ds + wr(t)* u1(t))/wr(t);
% ud(t) == (u2(t)*wr(t))/wd;


Eqns = [diff(h(t),t) == - u1(t)
        diff(u1(t),t) == g*(h(t)/ds) - br * (u1(t)); 
        diff(u2(t),t) == -g * (h(t)/ds) - br * (u2(t));
        diff(ud(t),t) == -g * (h(t)/ws) - bd * (ud(t));
        diff(wr(t),t) == -us];
        
[DEsys,Subs] = odeToVectorField(Eqns);
DEfcn = matlabFunction(DEsys, 'Vars',{T,Y});
tspan = linspace(0, 200, 251);
Y0 = [0; 0; 0; 0; 30]+0.001;
[T, Y] = ode45(DEfcn, tspan, Y0);
figure
for k = 1:size(Y,2)
    subplot(size(Y,2),1,k)
    plot(T, Y(:,k))
    grid
    legend(string(Subs(k)))
end

[Star Strider]

The differential equations (‘Eqns’) are the only ones used to create the vector field (‘DEsys’) and therefore ‘DEfcn’. The 3 equations as written are only defined as ‘u1(t)’, ‘u2(t)’ and ‘ud(t)’, not as the expressions. If you instead use single equal signs so they are assignments, the code fails because odeToVectorField cannot parse them.

If you want them included in ‘DEsys’, you must manually substitute them, so that ‘Eqns’ becomes:

Eqns = [diff(u1(t),t) == g*(h(t)/ds) - br * ((-us * ds + wr(t) * u2(t))/wr(t)); 
        diff(u2(t),t) == -g * (h(t)/ds) - br * ((us*ds + wr(t)* u1(t))/wr(t));
        diff(ud(t),t) == -g * (h(t)/ws) - bd * ((u2(t)*wr(t))/wd);
        diff(h(t),t) == - u1(t);
        diff(wr(t),t) == -us];

[DEsys,Subs] = odeToVectorField(Eqns);
DEfcn = matlabFunction(DEsys, 'Vars',{T,Y});

That will likely do what you want.

The complete, revised, code is now:

ds = 18.32;
ws = 45;
br = 0.0015;
bd = 0.0015;
rho = 1019;
g = 9.81;
us = 0.06;
massvessel = ds * ws * rho;
wd = 1.13;

syms u1(t) u2(t) ud(t) h(t) wr(t) T Y

Eqns = [diff(u1(t),t) == g*(h(t)/ds) - br * ((-us * ds + wr(t) * u2(t))/wr(t)); 
        diff(u2(t),t) == -g * (h(t)/ds) - br * ((us*ds + wr(t)* u1(t))/wr(t));
        diff(ud(t),t) == -g * (h(t)/ws) - bd * ((u2(t)*wr(t))/wd);
        diff(h(t),t) == - u1(t);
        diff(wr(t),t) == -us];

[DEsys,Subs] = odeToVectorField(Eqns);
DEfcn = matlabFunction(DEsys, 'Vars',{T,Y});

 tspan = linspace(0, 200, 251);
Y0 = [0; 0; 0; 0; 30]+0.001;
[T, Y] = ode45(DEfcn, tspan, Y0);

figure
for k = 1:size(Y,2)
    subplot(size(Y,2),1,k)
    plot(T, Y(:,k))
    grid
    title(string(Subs(k)))
end

The order would of course still be important.

.

[Iris Heemskerk]

Okay, that is what I thought the problem was. However if you plot this code and you look at the graph of ud(t).

The outcome of the solver (ud) does not agree with ud = wr * u2/wd (figure 2 in the code).

close all;
clear all;

ds = 18.32
ws = 45
br = 0.0015
bd = 0.0015
rho = 1019
g = 9.81
us = 0.06
massvessel = ds * ws * rho
wd = 1.13

syms u1(t) u2(t) ud(t) h(t) wr(t) T Y


Eqns = [diff(u1(t),t) == g*(h(t)/ds) - br * ((-us * ds + wr(t) * u2(t))/wr(t)); 
        diff(u2(t),t) == -g * (h(t)/ds) - br * ((us*ds + wr(t)* u1(t))/wr(t));
        diff(ud(t),t) == -g * (h(t)/ws) - bd * ((u2(t)*wr(t))/wd);
        diff(h(t),t) == - u1(t);
        diff(wr(t),t) == -us];
    
[DEsys,Subs] = odeToVectorField(Eqns);

DEfcn = matlabFunction(DEsys, 'Vars',{T,Y});

tspan = linspace(0, 400, 401);

Y0 = [0; 0.0675; 1.41; 0; 30]+0.001;

[T, Y] = ode45(DEfcn, tspan, Y0);

figure
for k = 1:size(Y,2)
    subplot(size(Y,2),1,k)
    plot(T, Y(:,k))
    grid
    legend(string(Subs(k)))
end

figure
plot(T, Y(:,3))
hold on
plot(T, Y(:,1).*Y(:,5)./wd)
legend('ud', 'ud = u2 * wr/wd')

[Star Strider]

It will not agree, because the value of that is plotted is the integrated value of as coded in ‘Eqns’.