How to convert a function handle to symbolic function?

9 vues (au cours des 30 derniers jours)
Vedant Kumar
Vedant Kumar le 1 Jan 2020
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);

Réponses (2)

Steven Lord
Steven Lord le 17 Sep 2024
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.

Pramil
Pramil le 17 Sep 2024
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.

Community Treasure Hunt

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

Start Hunting!

Translated by