problem in creating function file using ode45

I am trying to solve a second order coupled ode using ode45, however I am facing issues in making it function file and providing variables as input. How to convert symbolic input to numerical in main file using ode45.

Answer 1

[Torsten]

Be careful with the order of the solution vector Y. It is (y, Dy, x, Dx). You have to supply the vector of initial values "ic" in the same order !

rng("default")
clc
clear all
close all
O=rand;
a=rand;
g=9.81;
L=rand;

[y]=coupled_ode(O, a, g, L);

VF = 

Subs = 

ftotal = function_handle with value:
    @(t,Y)[Y(2);Y(1).*(-7.725211393067794e+1)-sin(9.057919370756192e-1).*Y(4).*1.629447372786358;Y(4);Y(3).*(-7.725211393067794e+1)+sin(9.057919370756192e-1).*Y(2).*1.629447372786358]

Disp=y;

function [y] = coupled_ode(Onum, anum, gnum, Lnum)
  syms O a g L x(t) y(t) t Y 
  dx = diff(x);
  d2x = diff(x,2);
  dy = diff(y);
  d2y = diff(y,2);
  Eq1 = d2x == 2*O*sin(a)*dy - (g/L)*x(t);
  Eq2 = d2y == -2*O*sin(a)*dx - (g/L)*y(t);
  [VF,Subs] = odeToVectorField(Eq1, Eq2)
  VF = subs(VF,[O a g L],[Onum anum gnum Lnum]);
  ftotal = matlabFunction(VF,'Vars',{t,Y})

  tspan = [0 25];                                             % Choose Appropriate Simulation Time
  ic = [0 1 0 1];                                             % Choose Appropriate Initial Conditions
  [t,y] = ode45(ftotal, tspan, ic);

  figure
  plot(t, y)
  grid
  legend(string(Subs))
end

[Susmita Panda]

Thank you @ Torsten sir. rng(default) helped me.This is exactly what I was searching.I didnt have any idea about this command.

[Torsten]

It's always a good idea if programming parts that serve different tasks have their own subfunctions. But in your case, I don't understand why and how you divided your code in the two different .m files.

E.g. to have control over what is done in the function coupled_ode, you will have to pass tspan, the initial conditions for the variables and the parameter values for O, a, g and L. The values returned should be t and y. The plotting part also should be done in the main program because it makes no sense to do the plotting in coupled_ode while passing the results as output parameters back to the calling program.

I'd arrange the code as

rng("default")
clc
clear all
close all

O = rand;
a = rand;
g = 9.81;
L = rand;

f = function_to_integrate(O, a, g, L)

f = function_handle with value:
    @(t,Y)[Y(2);Y(1).*(-7.725211393067794e+1)-sin(9.057919370756192e-1).*Y(4).*1.629447372786358;Y(4);Y(3).*(-7.725211393067794e+1)+sin(9.057919370756192e-1).*Y(2).*1.629447372786358]

tstart = 0.0;
tend = 25.0;
tspan = linspace(tstart,tend,100);                          % Choose Appropriate Simulation Time
ic = [0 1 0 1];                                             % Choose Appropriate Initial Conditions
[t,y] = ode45(f, tspan, ic);

figure
plot(t, y)
grid

legend(string(Subs))

Unrecognized function or variable 'Subs'.

function f = function_to_integrate(Onum, anum, gnum, Lnum)
  syms O a g L x(t) y(t) t Y 
  dx = diff(x);
  d2x = diff(x,2);
  dy = diff(y);
  d2y = diff(y,2);
  Eq1 = d2x == 2*O*sin(a)*dy - (g/L)*x(t);
  Eq2 = d2y == -2*O*sin(a)*dx - (g/L)*y(t);
  [VF,Subs] = odeToVectorField(Eq1, Eq2);
  VF = subs(VF,[O a g L],[Onum anum gnum Lnum]);
  f = matlabFunction(VF,'Vars',{t,Y});
end

Answer 2

[Alan Stevens]

Are you looking for something like this? (I might have misinterpreted your equations!)

O=1;
a=pi/2;
g=9.81;
L=100;

tspan = [0 15];                                             % Choose Appropriate Simulation Time
ic = [0 1 0 1];                                             % Choose Appropriate Initial Conditions
[t,y] = ode45(@(t,y) ftotal(t,y,O,a,g,L), tspan, ic);

figure
plot(t, y)
grid
legend('y1','dy1/dt','y2','dy2/dt')

function dydt=ftotal(~,y,O, a, g, L)

y1 = y(1); v1 = y(2); y2 = y(3); v2 = y(4);

dydt = [ v1;
         2*O*sin(a)*y2 - g/L;
         v2;
         -2*O*sin(a) - g/L*y1];
end

[Susmita Panda]

@Alan sir, this is a demo question, I have two ode equations coupled in such a such that the acceleration of one variable x also contains y also, e.g xddot=yddot+A+B ydot= xdot+M+N

I have solved successfully using the procedure however i am facing issue in providing variables like in this question.

Is there any way we can defined acceleration of first variable with the second variable using your procedure?

[Alan Stevens]

Hmm! You need an equation for yddot. What are your actual equations? Try uploading them exactly in mathematical notation (using an image if necessary).

[Susmita Panda]

Sir please find the attached file containing second order coupled odes. Is there any simple way to do using ode45 because the method I tried having problem of defining variables in another function file.

[Torsten]

Call ode45 with the function

 fun = @(t,y) [y(2); xddot; y(4); yddot]

where xddot and yddot are computed as

syms a b c d e f g h i t x y xdot ydot xddot yddot 
eqn1 = (1-e)*xddot + e*yddot == -2*a*xdot - x - sin(b*t) + c*(x-y) + d*(xdot-ydot)

eqn1 = 

eqn2 = (g-h)*xddot - (i+g)*yddot == -f*(x-y)

eqn2 = 

solve([eqn1 eqn2],[xddot yddot])

ans = struct with fields:
    xddot: -(g*x + i*x + g*sin(b*t) + i*sin(b*t) + 2*a*g*xdot - c*g*x + 2*a*i*xdot + c*g*y + e*f*x - d*g*xdot - c*i*x - e*f*y + d*g*ydot + c*i*y - d*i*xdot + d*i*ydot)/(g + i - e*h - e*i)
    yddot: (f*x - f*y - g*x + h*x - g*sin(b*t) + h*sin(b*t) - 2*a*g*xdot + 2*a*h*xdot + c*g*x - c*g*y - c*h*x - e*f*x + d*g*xdot + c*h*y + e*f*y - d*g*ydot - d*h*xdot + d*h*ydot)/(g + i - e*h - e*i)

[Susmita Panda]

Thank you @Torsten sir for this method. I will try this.

but is there any way we can call function "coupled_ode" in "Main_file" which I posted without error. This method is convenient to me and produces good result when they variables are called in same file.I am getting error of this kind when I am making it function file and providing variables in "Main_file":

Error using superiorfloat

Inputs must be floats, namely single or double.

Error in odearguments (line 114)

dataType = superiorfloat(t0,y0,f0);

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in coupled_ode (line 18)

[t,y] = ode45(@(t,y) ftotal(t,y,O,a,g,L), tspan, ic);

Error in Main_file (line 9)

[y]=coupled_ode(O, a, g, L);

[Walter Roberson]

After you get the xddot and yddot using solve(), call odeFunction to create something to pass to ode45() .

I recommend following the workflow shown in the first example to odeFunction