Solving set of differential equations

I made these equations on simulink and their outputs are correct, but making them on .m file is really mess for me

y,z,u are my dependent variables.

x is 0.007*sin(4*pi*t)

so i solve EQ5 then subs in 3 after that i solve 2&4 together and subs back into 1 to get f

but plotting f with time is really wrong... this is the code

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tspan = [0 2];

u0 = 0; v=1.5; eta=190; n=2 ; co_a=784; co_b=1803; ko=3610;

c1_a=14649; c1_b=34622; k1=840; xo=0.0248; alpha_a=1241; alpha_b=38430;

beta=2059020; gamma=136320; A=58;

%%%%%%% Eqn 2 & 4

yZero = [0,0]; % and an initial condition

[T,myY]=ode45(@Eqn2_4,tspan,yZero); % solve

[t,u] = ode45(@(t,u) -eta*(u-v), tspan, 0);

y=myY(:,1);

z=myY(:,2);

k1=840;

x= 0.007*sin(4*pi*t);

xdot= 0.007*4*pi*cos(4*pi*t);

y = interp1(T,y,t);

z = interp1(T,z,t);

c1=c1_a+c1_b.*u; %%eqn3

co=co_a+co_b.*u; %%eqn3

alpha=alpha_a+alpha_b.*u; %%eqn3

ydot=(1./(co+c1).*((alpha.*z)+(co.*xdot)+ko.*(x-y)));

fa=c1.*(ydot)+k1*(x-xo);

plot(t,fa)

%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

and the function Eqn2_4 is

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ret]=Eqn2_4(t,y)

% set up ret in correct form (i.e. a column vector)

ret=zeros(2,1);

u0 = 0; tspan = [0 2];

v=1.5;

eta=190; n=2 ;

co_a=784; co_b=1803; ko=3610;

c1_a=14649; c1_b=34622; k1=840;

xo=0.0248; alpha_a=1241; alpha_b=38430;

beta=2059020; gamma=136320; A=58;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[t,u] = ode45(@(t,u) -eta*(u-v), tspan, u0);

alpha=alpha_a+alpha_b.*u; %%eqn3

c1=c1_a+c1_b.*u; %%eqn3

co=co_a+co_b.*u; %%eqn3

% extract current values of V and P

y1 = y(1);

z = y(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%

alpha=interp1(u,alpha,t);% Interpolate the data set

c1 = interp1(u,c1,t); % Interpolate the data set (u,c1) at time t

co = interp1(u,co,t); % Interpolate the data set (u,co) at time t

x= 0.007*sin(4*pi*t);

xdot= 0.007*4*pi*cos(4*pi*t);

% return updates

ret(1) = (1./(co(1)+c1(1)).*((alpha(1).*z)+(co(1).*xdot(1))+ko.*(x(1)-y1)));

ret(2) = (-gamma.*abs(xdot(1)-ret(1)).*abs(z))-(beta.*(xdot(1)-ret(1)).*pow2(abs(z)))+50*(xdot(1)-ret(1));

end

[Torsten]

Solve (2),(4) and (5) together in one call to ode15i.

[Mimo T]

Thanks i am trying it

Answer 1

[William Rose]

@mohamed talha,

@Torsten is exactly right, and @Torsten is incredibly good at quickly seeing the key thing to do.

In case it is not obvious, you must re-write your equations using two vectors, let's call them y and yp (for y prime)

,

where y and yp on the left are the new y and yp, which you will pass to ode15i, and the y,z,u on the right are the old names for the variables you want to solve for.

Rewrite equations 2, 4, 5, using the new y and yp. As you do this, replace x with 0.007*sin(4*pi*t) and replace xdot with 0.028*pi*cos(4*pi*t), and replace α with , and replace with , and replace with , and write equations 2, 4, 5 so that you have zero on one side and everything else on the other side. Then equations 2, 4, 5 become

(You can finish equation 4.)

Now the equations are in the form that you need in order to write the function you will pass to ode15i:

function res = mydiffeq(t,y,yp)
%MYDIFFEQ returns residual of my implicit differential equations
eta=190; n=2; k0=3610;
c0a=784; c0b=1803; 
c1a=14649; c1b=34622; 
alphaa=1241; alphab=38430;
beta=2059020; gamma=136320;
res=[yp(1)-((alphaa+alphab*y(3))*y(2)+(c0a+c0b*y(3))*0.028*pi*cos(4*pi*t)+k0*(0.007*sin(4*pi*t)-y(1)))/(c0a+c0b*y(3)+c1a+c1b*y(3));...
    yp(2)+gamma*abs(0.028*pi*cos(4*pi*t)-yp(2))*abs(y(2))^(n-1)*y(2)+beta*(...);...
    yp(3)+eta*(y(3)-nu));
end

Read the help for ode15i, especially the examples for the Weissinger problem and the Robertson problem. THey will show you how to use ode15i.

Check all my equations and code above for possible errors. Try it. Good luck.

[Mimo T]

@William Rose thanks alot, but i don't get what to write in y and z vectors ?!

[Torsten]

Your solution variables in William's code are put in arrays called y and yp.

They correspond with your y, z and u variables as follows:

y(1) is your y
yp(1) is your ydot
y(2) is your z
yp(2) is your zdot
y(3) is your u
yp(3) is your udot

Now compare and try to recover your 3 equations from William's "res" vector:

res = [yp(1)-((alphaa+alphab*y(3))*y(2)+(c0a+c0b*y(3))*0.028*pi*cos(4*pi*t)+k0*(0.007*sin(4*pi*t)-y(1)))/(c0a+c0b*y(3)+c1a+c1b*y(3));...
       yp(2)+gamma*abs(0.028*pi*cos(4*pi*t)-yp(2))*abs(y(2))^(n-1)*y(2)+beta*(...);...
       yp(3)+eta*(y(3)-nu));

[William Rose]

@Torsten, thank you.

[William Rose]

@mohamed talha, I notice an error in the code I posted. The penultimate line of the function is missing a square bracket which is required to complete the vector res. Therefore the line should be

yp(3)+eta*(y(3)-nu))];

I try to be correct, but there could be other errors, so check my equations and code carefully.

[Mimo T]

@Torsten@William Rose first,thanks for this the answer that worked, for my first eqn how can i git my y' to get F.. should i differentiate my y data from ode15i or it is better to subs back in eqn 2 to get y' and use it in eqn 1 ?

[William Rose]

@mohamed talha, you're welcome. You have : it is called yp(:,1) in the new notation. It has been computed at the times in vector t, which you got back from ode15i. And c1=c1a+c1b*u. And "u" is called y(:,3) in your code. Therefore you can compute

f_a=(c1a+c1b*y(:,3)).*yp(:,1)+k1*(0.007*sin(4*pi*t)-x0);

and it will be valued at the times in t.

[Mimo T]

@William Rose

[t,myY]=ode15i(@mydiffeq,[0 0.1],[0; 0; 0],[0;0;0]);

i don't have access to my yp ! i only get three colums output from ode15i which are y z u

[William Rose]

You are corrct. My mistake. I would compute numerically:

%Next: Compute a first order estimate of the derivative of y(1,:)
ydot=diff(y(:,1))./diff(t);  
%Next: compute f_a
f_a=(c1a+c1b*y(2:end,3)).*ydot+k1*(0.007*sin(4*pi*t(2:end))-x0);

I use y(2:end,3) and t(2:end) because diff(y(:,1)) and diff(t) returns N-1 elements, if y and t have N elements. Therefore ydot will have N-1 elements. Try it.

[Torsten]

@mohamed talha

After integration, use "deval" to get the derivatives of y, z and u with respect to time:

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

Here is a similar question:

https://de.mathworks.com/matlabcentral/answers/216540-what-is-the-most-elegant-way-to-get-yp-from-solution-of-ode15i