preparing symbolic ODE for ode45

Consider the following symbolic ODE:

syms gamma(t)
C1 = 2.4
C2 = 3.1
ode = diff(gamma, t, 2) == C1 * sin(gamma) + C2 * cos(diff(gamma))

Now, I want to convert the above expression to the function odefun which I can use with ode45:

function res = odefun(t, gamma)
    res = [gamma(2);
           2.4*sin(gamma(1)) + 3.1*cos(gamma(2))]
end

Are there ways to automate (parts of) this process? The above ode is only a simple example to illustrate my problem. The ode I'm really trying to solve has way more terms with coefficients.

Thanks

Answer 1

[Star Strider]

Try something like this —

syms gamma(t) T Y 
C1 = 2.4

C1 = 2.4000

C2 = 3.1

C2 = 3.1000

ode = diff(gamma, t, 2) == C1 * sin(gamma) + C2 * cos(diff(gamma))

ode(t) = 

[VF,Subs] = odeToVectorField(ode)

VF = 

Subs = 

odefun = matlabFunction(VF, 'Vars',{T,Y})

odefun = function_handle with value:
    @(T,Y)[Y(2);cos(Y(2)).*(3.1e+1./1.0e+1)+sin(Y(1)).*(1.2e+1./5.0)]

ic = [0 1];
tspan = [0 10];

[t,omega] = ode45(odefun, tspan, ic);

figure
plot(t, omega)
grid
legend(string(Subs), 'Location','best')

The constants do not have to be specified in the original equations. They can be included as parameters, so for example —

syms gamma(t) T Y C1 C2

ode = diff(gamma, t, 2) == C1 * sin(gamma) + C2 * cos(diff(gamma))

odefun = matlabFunction(VF, 'Vars',{T,Y, [C1,C2]})

ic = [0 1];

tspan = [0 10];

CV = rand(1,2)

[t,omega] = ode45(@(t,omega)odefun(t,omega,CV), tspan, ic);

Now they can be changed in the numeric code (for example in a loop or nested loops) without having to re-derive them each time in the original symbolic code.

.

[Clemens Herrmann]

Thanks a lot, works like a charm!

[Star Strider]

As always, my pleasure!

.