**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# How to solve two coupled differential equations using ode45.

92 views (last 30 days)

Show older comments

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.

##### 0 Comments

### Accepted Answer

Alan Stevens
on 19 Oct 2020

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

##### 22 Comments

Musanna Galib
on 19 Oct 2020

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
on 19 Oct 2020

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
on 19 Oct 2020

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
on 19 Oct 2020

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
on 19 Oct 2020

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?

Musanna Galib
on 20 Oct 2020

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
on 20 Oct 2020

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
on 20 Oct 2020

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

Musanna Galib
on 21 Oct 2020

Alan Stevens
on 21 Oct 2020

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
on 21 Oct 2020

Alan Stevens
on 21 Oct 2020

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
on 21 Oct 2020

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
on 21 Oct 2020

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
on 21 Oct 2020

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
on 22 Oct 2020

haohaoxuexi1
on 19 Jan 2022

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
on 19 Jan 2022

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
on 19 Jan 2022

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.

### More Answers (1)

Alan Stevens
on 20 Jan 2022

"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.

##### 0 Comments

### See Also

### Tags

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.