How to convert a function handle to symbolic function?

Hello!

I first performed curve fit on a dataset with fourier8 option and then converting the cfit object into a function handle. I intend to use this function handle in dsolve but couldn't figure out how to convert a function handle to symbolic function. The aim to use dsolve is to obtain exact solution of an ODE, so any other alternative method which does the trick will also be very appreciated. Attaching the code and dataset for your kind reference:

clear all
load 'JFM_angular.mat'
fit_JFM=fit(time_JFM,C_moment_JFM,'fourier8');
fun_JFM=@(t) fit_JFM(t);%creating a function for the moment coefficient
fun_JFM_sym = sym(fit_JFM);

I=4.1;%moment of inertia
U_r=6.6;%reduced velocity
k=4*pi^2*I/U_r^2;%spring stiffness
damping=0.15;% percentage of critical damping
b=2*damping*sqrt(k*I);%damping parameter
theta_o=0;%initial angle
syms theta(t)
Dtheta=diff(theta,t);%first derivative of theta
D2theta=diff(theta,t,2);%second derivative of theta
ode= I*D2theta + b*Dtheta + k*(theta-theta_o) == fun_JFM_sym; %writing the ODe

cond1 = theta(0) == theta_o;
cond2 = Dtheta(0) == 0;%zero velocity at time=0
cond=[cond1 cond2];
theta_sol(t)=dsolve(ode,cond);

Answer 1

[Steven Lord]

Are you hoping to use this fit object as the right-hand side of a differential equation ? If so keep in mind that fit objects have a number of methods that you can call on them, one of which is integrate.

x = 0:5;
y = 2*x.^2 + 3*x - 4;
yint = (2/3)*x.^3 + (3/2)*x.^2 - 4*x; % integral of y, computed manually

x2 = 0:0.25:5;
f = fit(x.', y.', 'poly2'); % create a sample fit object
integralOverTime = integrate(f, x2, 0);

plot(x2, f(x2), '-o', DisplayName = "fitted curve");
hold on
plot(x, y, '+', DisplayName = "original points")
legend show

figure
plot(x2, integralOverTime, '-o', DisplayName = "integral of fit")
hold on 
plot(x, yint, '+', DisplayName = "manually computed integral")
legend show

Those curves look to be in pretty good agreement.

Answer 2

[Pramil]

Hey Vedant,

To resolve the issue, you need to convert the “cfit” object from the curve fitting process into a symbolic expression that can be used in “dsolve”.

Using “sym(fit_JFM)” directly won't work because “fit_JFM” is not directly convertible to a symbolic expression. Instead, you can manually extract the coefficients from the Fourier fit and construct a symbolic expression.

Here's how you can achieve this:

1. Extract the coefficients from the “cfit” object.
2. Construct the symbolic expression using these coefficients.
3. Use this symbolic expression in your differential equation.

Here’s how the code would be, after applying the above changes:

clear
load 'JFM_angular.mat'
fit_JFM = fit(time_JFM, C_moment_JFM, 'fourier8');

% Extract coefficients from the fit
coeffs = coeffvalues(fit_JFM);
w = coeffs(1); % fundamental frequency

% Construct symbolic Fourier series expression
syms t
fourier_expr = coeffs(2);
for n = 1:8
    fourier_expr = fourier_expr + coeffs(2*n+1) * cos(n*w*t) + coeffs(2*n+2) * sin(n*w*t);
end

I = 4.1; % moment of inertia
U_r = 6.6; % reduced velocity
k = 4*pi^2*I/U_r^2; % spring stiffness
damping = 0.15; % percentage of critical damping
b = 2*damping*sqrt(k*I); % damping parameter
theta_o = 0; % initial angle

syms theta(t)
Dtheta = diff(theta, t); % first derivative of theta
D2theta = diff(theta, t, 2); % second derivative of theta

% Set up the ODE using the symbolic Fourier expression
ode = I*D2theta + b*Dtheta + k*(theta - theta_o) == fourier_expr;

cond1 = theta(0) == theta_o;
cond2 = Dtheta(0) == 0; % zero velocity at time=0
conds = [cond1 cond2];

theta_sol(t) = dsolve(ode, conds);

disp(theta_sol)

Hope it helps.