How to solve two coupled differential equations using ode45.

I need help to solve two paired differential equations. After hours of unsuccessful attemps, I am asking for help from you. Any suggestion will be highly appriciated.

The equations are -

m*x''(t) + U*x(t) + γ*x'(t) + kv*V(t) - ζ = 0

V'(t) - kc*x'(t) + τ*V(t) = 0

initial conditions are x(0) = 1, X'(0) = 0, V(0) =0

x and V are the variables. All others are constant.

Thanks in advance.

Answer 1

[Alan Stevens]

The following should help:

% tspan = [0 tend];
% IC = [1 0 0];      % initial conditions 
% [t, X] = ode45(@fn, tspan, IC);
% x = X(:,1);
% z = X(:,2);
% V = X(:,3);
%
% function dXdt = fn(t,X)
%           x = X(1); z = X(2); V = X(3);
%           constants = ... list them with their values
%           dXdt = [z;
%                   -(y*z+U*x+kv*V)/m;
%                   kc*z-tau*V];
% end        

[Musanna Galib]

Thanks for your reply. I successfully solved this. Thanks a lot.

Can you kindly help in another problem? I have to solve the problem for different values of a constant Delta. (for example Delta = 0.1: 0.1: 10) . But I am getting this error for elementwise calculation in ode45. --"Error using vertcat - Dimenions of arrays being concatenated are not consistent "

Delta = 0.1:0.1:1;

f = @(t,x) [x(2); (1/m)*((-2*k*x(1)-1.5*(A*x(1)^2+B.*Delta.^2).^-2.5) ...
    - gamma*x(2) - K_v*x(3)); K_c*x(2)-((1/tau_p)*x(3))]  
[ts,xs] = ode45(@(t,x)f(t,x),[t0,T],initial_conditions);    

[Alan Stevens]

Like this?

% Replace the following constants with your own values
 m = 1; k = 1; A = 1; B = 1; gamma = 1; K_v = 1; K_c = 1; tau_p = 1;
 t0 = 0; T = 0.5;
 initial_conditions = [1 0 0];
 % Define function f before defining Delta or Matlab thinks you want all
 % the Delta values in the function at the same time. Make f a function of
 % Delta.
 
f = @(t,x,Delta) [x(2); (1/m)*((-2*k*x(1)-1.5*(A*x(1)^2+B.*Delta.^2).^-2.5) ...
    - gamma*x(2) - K_v*x(3)); K_c*x(2)-((1/tau_p)*x(3))];

Delta = 0.1:0.1:1;
% Loop through the values of delta
for i = 1:numel(Delta)
[ts, xs] = ode45(@(t,x)f(t,x,Delta(i)),[t0:T/100:T],initial_conditions);  
xs1(:,i) = xs(:,1);
xs2(:,i) = xs(:,2);
xs3(:,i) = xs(:,3);
end
subplot(3,1,1)
plot(ts,xs1),grid
subplot(3,1,2)
plot(ts,xs2),grid
subplot(3,1,3)
plot(ts,xs3),grid

[Musanna Galib]

Thank you very much Alan. It worked well.

Another question if you don't mind. I have to draw P with respect to Delta. Can't figure it out how to relate P with delta.

Here,

P = xs(:,3).^2./R_L;

[Alan Stevens]

Stick this on the end of the previous coding

P = xs3.^2/R_L;
% If you want to plot P vs Delta for a specific times, 
% say t = T/2 and t = T then
figure(2);
plot(Delta, P(51,:),'r',Delta,P(end,:),'b'),grid
xlabel('Delta'),ylabel('P')
legend(['time = ', num2str(ts(51))],['time = ',num2str(ts(end))])

% If you want a plot for all times and all Delta's then perhaps
% a 3D surf plot is better
[Dx, Dy] = meshgrid(Delta,ts);
figure(3)
surf(Dx,Dy,P)
xlabel('Delta'),ylabel('ts'),zlabel('P')

[Musanna Galib]

Thank you very much for the reply. I think some modification is needed in the code. P should have same number of columns as xs3. But it is not the case here.

Also the graph represents a straight line even though it should be a curve according to the data!

Can you kindly check?

[Alan Stevens]

When I ran the code P had the same number of columns as xs3.

[Musanna Galib]

Sorry for the late reply. Yes, you are correct. The code is ok. Thanks. I would like to ask for another help.

I want to add a white gaussian noise in my system with mean = 0, variance, σ= 3,6,12*10^-4 and exponential autocorrelation function with correlation time 𝜏 =0.1. Can you help me with this?

[Alan Stevens]

I'm no expert on those aspects, but for the gaussian noise type doc randn in the command window.

For exponential autocorrelation have a look at: https://uk.mathworks.com/help/signal/ug/autocorrelation-function-of-exponential-sequence.html

[Musanna Galib]

Thank you very much. I will look into this.

[Musanna Galib]

Hello Alan, I am facing an error in my code when I changes the values of Keff to 7.11*10^-3 from 1. Can you kindly look into it?

% Constants
K_c = 0.5;
omega = 10;
gamma = 0.016; 
K_v =1;
eta_z = random('Norm',0,3.4*10^-4);
m = 0.0155;  
keff = 7.11*10^-3; 
d = 2.97;
mu_o = 12.57*10^-7; 
M = 0.051; 
R_L = 100*10^6; 
C = 112*10-9; 
tau_p = R_L*C;

%functions definition

k = keff/2;
A = d^2*(mu_o*M^2/(2*pi*d))^(-2/3);
B = A/d^2;
C = k/d^2;
U = k*x(1)^2 + ((A*x(1)^2)+(B*Delta.^2)).^(-3/2) + C*Delta.^2
V = tau_p*(K_c*x(2) - V(1))


f = @(t,x,Delta) [x(2); (1/m)*((-2*k*x(1)-1.5*(A*x(1)^2+B.*Delta.^2).^-2.5) ...
    - gamma*x(2) - K_v*x(3) + 0.083*sin(omega*t)); ...
     K_c*x(2)-((1/tau_p)*x(3))]    


t0 = 0;                            
initial_conditions = [0;0;0];      
T = 100;                           
Delta = .005:0.001:0.01;

 % solve IVP
for i = 1:numel(Delta)
[ts,xs] = ode45(@(t,x)f(t,x,Delta(i)),[t0,T],initial_conditions);    
xs1(:,i) = xs(:,1);
xs2(:,i) = xs(:,2);
xs3(:,i) = xs(:,3);
end



fprintf('x(T) = %g, x''(T) = %g\n',xs(end,1),xs(end,2))

disp('       t           x1          x2         x3=V')
disp([ts,xs])                     

plot(ts,xs(:,1),'b')              
title('Solution x(t) of IVP')
xlabel('t'); 
ylabel('x(t)');grid on


plot(ts,xs(:,3))                  %plot V
title('Solution V(t) of IVP')
xlabel('t'); 
ylabel('V(t)');grid on


U = k.*xs(:,1).^2 + (A.*xs(:,1).^2+(B.*Delta.^2)).^(-3/2) + C.*Delta.^2   
plot(xs(:,1),U(:,:),'b')          
title('Solution U(t) of IVP')
xlabel('x'); 
ylabel('U(x)');grid on


P = xs3.^2/R_L;
plot(Delta, P(50,:),'r',Delta,P(end,:),'b'),grid
xlabel('Delta'),ylabel('P')
legend(['time = ', num2str(ts(50))],['time = ',num2str(ts(end))])
plot(ts,P(:,:),'b')           
title('Electric Power as a function of time')
xlabel('t'); 
ylabel('P');grid on

[Alan Stevens]

Will do, but it will have to wait until tomorrow now!

[Musanna Galib]

No problem at all!. The problem I am facing is whenever I change values of variables it sayes,"Unable to perform assignment because the size of the left side is 3829 by 1 and the size of the right side is 3825 by 1"

[Alan Stevens]

The problem arises because, left to its own devices ode45 selects and returns results at times of its own choosing. This means that for one value of "i" in the loop there could be a different number of results returned from that of another. Therefore, you need to force ode45 to return the same number of resuts every sweep through the loop. Thats done in the following (I've also made a couple of other notes in the listing):

% Constants
K_c = 0.5;
omega = 10;
gamma = 0.016; 
K_v =1;

% eta_z = random('Norm',0,3.4*10^-4);
% I don't have the Statistics and Machine Learning toolbox, so
% I've replaced the above with:
eta_z = 3.4*10^-4*randn; 

m = 0.0155;  
keff = 7.11*10^-3; 
d = 2.97;
mu_o = 12.57*10^-7; 
M = 0.051; 
R_L = 100*10^6; 
C = 112*10-9; 
tau_p = R_L*C;

%functions definition

k = keff/2;
A = d^2*(mu_o*M^2/(2*pi*d))^(-2/3);
B = A/d^2;
C = k/d^2;

% x isn't defined yet so I've commented out the following two lines
%U = k*x(1)^2 + ((A*x(1)^2)+(B*Delta.^2)).^(-3/2) + C*Delta.^2;  
%V = tau_p*(K_c*x(2) - V(1));


f = @(t,x,Delta) [x(2); (1/m)*((-2*k*x(1)-1.5*(A*x(1)^2+B.*Delta.^2).^-2.5) ...
    - gamma*x(2) - K_v*x(3) + 0.083*sin(omega*t)); ...
     K_c*x(2)-((1/tau_p)*x(3))];  


t0 = 0;                            
initial_conditions = [0;0;0];      
T = 100;                           
Delta = .005:0.001:0.01;

% Left to its own devices ode45 selects and returns results at times of its
% own choosing.  This means that for one value of "i" below there 
% could be a different number of results returned from that of another.
% We can force ode45 to return the same number of resuts for every "i" 
% by specifying them as follows:
tspan = 0:T/100:T;       
 % solve IVP
for i = 1:numel(Delta)
[ts,xs] = ode45(@(t,x)f(t,x,Delta(i)),tspan,initial_conditions);    
xs1(:,i) = xs(:,1);
xs2(:,i) = xs(:,2);
xs3(:,i) = xs(:,3);
end



fprintf('x(T) = %g, x''(T) = %g\n',xs(end,1),xs(end,2))

disp('       t           x1          x2         x3=V')
disp([ts,xs])
% Use the figure command to separate the plots, otherwise 
% successive plots will replace each other.
figure(1)
plot(ts,xs(:,1),'b')              
title('Solution x(t) of IVP')
xlabel('t'); 
ylabel('x(t)');grid on

figure(2)
plot(ts,xs(:,3))                  %plot V
title('Solution V(t) of IVP')
xlabel('t'); 
ylabel('V(t)');grid on

figure(3)
U = k.*xs(:,1).^2 + (A.*xs(:,1).^2+(B.*Delta.^2)).^(-3/2) + C.*Delta.^2 ;  
plot(xs(:,1),U(:,:),'b')          
title('Solution U(t) of IVP')
xlabel('x'); 
ylabel('U(x)');grid on

figure(4)
P = xs3.^2/R_L;
plot(Delta, P(50,:),'r',Delta,P(end,:),'b'),grid
xlabel('Delta'),ylabel('P')
legend(['time = ', num2str(ts(50))],['time = ',num2str(ts(end))])
plot(ts,P(:,:),'b')           
title('Electric Power as a function of time')
xlabel('t'); 
ylabel('P');grid on

I have no idea if the results are sensible or not - that's for you to decide!

I should add that I set the time intervals to be returned by ode45 at multiples of T/100 for testing purposes. However, you might want to set that at something like T/5000 or so in order to pick up the higher frequency components.

[Musanna Galib]

Thank you very much for your kind help. I think the graph of U is not correct. The corresponding values of x and U should be arranged from assending order before plot. Can you help with that?

[Alan Stevens]

Under figure(3) do the following

U = k.*xs1.^2 + (A.*xs1.^2+(B.*Delta.^2)).^(-3/2) + C.*Delta.^2 ; 
[xx, id] = sort(xs(:,1));
plot(xx,U(id),'b')          

This results in

Looking more carefully at the results I can see that the values of xs (as stored in xs1, xs2 and xs3) do not seem to change with the different values of Delta, at least in the range you've chosen!!

[Musanna Galib]

Thank you very much. U seems ok now. Yes you are correct. The range of Delta should be different. I am feeling the same problem for P also. It should not be straight line but should be a curve. But couldn't understand how to sort it (maybe in assending order?). Didn't understand how to link P with Delta in this perspective.

Also can you help how to create automatic legend for the delta range. I mean for delta=0.05:0.01:0.1, there are 5 curves for U. How to give legend to them?

[Alan Stevens]

I overlooked the fact that there should be a figure 5:

figure(3)
U = k.*xs1.^2 + (A.*xs1.^2+(B.*Delta.^2)).^(-3/2) + C.*Delta.^2 ; 
[xx, id] = sort(xs(:,1));
plot(xx,U(id),'b')          
title('Solution U(t) of IVP')
xlabel('x'); 
ylabel('U(x)');grid on

figure(4)
P = xs3.^2/R_L;
plot(Delta, P(50,:),'r',Delta,P(end,:),'b'),grid
xlabel('Delta'),ylabel('P')
legend(['time = ', num2str(ts(50))],['time = ',num2str(ts(end))])

figure(5)
plot(ts,P(id),'b')           
title('Electric Power as a function of time')
xlabel('t'); 
ylabel('P');grid on

There is no point in plotting anything against Delta until you use a range such that xs changes with Delta. Figure 4 shows that the two values you have plotted don't change as Delta changes.

[Musanna Galib]

I actaully found the way to plot P. It will be the avarage of each column of P against Delta. Can you help with that?

Also can you help how to create automatic legend for the delta range. I mean for delta=0.05:0.01:0.1, there are 5 curves for U. Is there any way other then manually put them?

[Alan Stevens]

To get the average of each column of P simple use mean(P). At present there is no noticeable variation with Delta. Find a set of Deltas for which P changes or there is no pint in having a legend as you will only see one horizontal line (they will all plot on top of each other!).

[haohaoxuexi1]

@Alan Stevens

Hi Alan, I saw you guys conservation above, it really helps me a lot to understand this problem.

I have an equation like below.

Apart from y1 y2 Vr, the rest of them are constant.

I wrote my function like this

with a mass matrix

M=diag([1 0 1 1 1]);

dx(1)=x(2);
dx(2)=k1*(x(3)-x(1))-(knl_1*x(1)+knl_3*x(1)^3);
dx(3)=x(4);
dx(4)=(-c*x(4)-k1*(x(3)-x(1))+kAmp*cos(w*t)+theta_p*x(5))/m;
dx(5)=(-theta_p*(kspring/(kpiezo+kspring))*(x(4)-x(2))*R_s-x(5))/R_s/C_p;

The reason why the second term in the mass matrix is zero, is because the dx(2) is a algebra equation. I am using ode23t to solve the problem.

But it also gives me an error like "This DAE appears to be of index greater than 1."

Would you mind to give me some suggestion how to solve it?

Thanks.

[Alan Stevens]

I'd be inclined to manipulate the equations to get the following (I'm assuming you know time zero values for y1, y2, dy2/dt and V2):

% Let dy2/dt = v2

% Define the following functions
% fna = @(t,y1,y2,v2,VR) (F0*cos(omega*t)+theta*VR+keq*(y1-y2)-c*v2)/m;
% fnb = @(y1,v2) keq*v2/(knl1+3*knl2*y1^2+keq);


% Set up a rate of change function that will be called by an ODE solver
% function dYdt = rate(t,Y,fna,fnb,fnc)
%          y1 = Y(1);
%          y2 = Y(2);
%          v2 = Y(3);
%          VR = Y(4);
%          dYdt = [fnb(y1,v2);
%                  v2;
%                  fna(t,y1,y2,v2,VR);
%                  (-R*theta*k1/(k1+kp)*(v2-fnb(y1,v2))-VR)/(R*Cp)];
% end

I can't check these numerically as you haven't supplied vaies for the constants, initial conditions or integration time.

If you intend to try this approach you must very carefully double check that I have manipulated the equations correctly!

[haohaoxuexi1]

dx(1)=(k1*x(3))/(knl_1+3*knl_3*x(1)^2+k1);
dx(2)=x(3);
dx(3)=(-c*x(3)-k1*(x(2)-x(1))+kAmp*cos(w*t)+theta_p*x(4))/m;
dx(4)=(-theta_p*(kspring/(kpiezo+kspring))*(x(3)-((k1*x(3))/(knl_1+3*knl_3*x(1)^2+k1)))*R_s-x(4))/R_s/C_p;

I changed the code to the above way and removed the mass matrix anymore, which follow your suggestions, and it worked with ode45.

I want to know if there is any logical mistake in my previous code? Why is my previous code is not working properly by treating the problem to be a DAE.

Also I noticed you wrote the fna fnb as separate function in your code instead of putting them all in dYdt function? Is there specific reason for you to do it this way? Like it can help calculating faster or it is just your habit?

Thank you for your help anyway.

Answer 2

[Alan Stevens]

"Why is my previous code is not working properly by treating the problem to be a DAE."

I don't know!

"Also I noticed you wrote the fna fnb as separate function in your code instead of putting them all in dYdt function? Is there specific reason for you to do it this way? Like it can help calculating faster or it is just your habit?"

You could put them inside the dYdt function here if you like. I tend to define them outside in case I need to call them separately outside of dYdt.

[haohaoxuexi1]

@Alan Stevens

Is there any way to solve it analytically? Or any sample code can be shared