How to convert a function handle to symbolic function?

10 次查看（过去 30 天）
Vedant Kumar
Vedant Kumar 2020-1-1
回答： Steven Lord 2024-9-17
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);

请先登录，再回答此问题。

回答（2 个）

Steven Lord
Steven Lord 2024-9-17
Ran in:
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 2024-9-17
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.

类别

在 Help Center 和 File Exchange 中查找有关 Equation Solving 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by