You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How to approximate a multivariable arctan function?
13 views (last 30 days)
Show older comments
Hi,
let the following function be given:
Is there a way in Matlab to approximate this function as a multivariable (x1,x2,x3,u1) rational polynomial?
lv and bv are positive constants.
x2 should be defined within -pi/4 and pi/4 [rad].
x3 should be defined within (-pi/4) and (pi/4) [rad/sec].
x1 > 0 [meter/sec].
u1 should be defined within (-pi/12) and (pi/12) [rad].
Thanks a lot in advance!
19 Comments
Jannis Erz
on 18 Sep 2022
I'd like to represent the function as a sum-of-squares to perform lyapunov stability analysis.
I need the function to be represented at least as a rational polynomial including four variables to be able to do that.
Jannis Erz
on 18 Sep 2022
This was just an example of a rational model! I'm looking for a rational model incl. 4 independent variables of course!
Torsten
on 18 Sep 2022
Edited: Torsten
on 18 Sep 2022
Ok.
Then evaluate your function for a number of points over the relevant intervals and make a least-square fit for the coefficient of the rational function you want to use.
A MATLAB program you could use for this task is lsqcurvefit.
But atan has a very restricted codomain. I wonder how rational functions in 4 variables should be able to approximate it.
Sam Chak
on 18 Sep 2022
Edited: Sam Chak
on 18 Sep 2022
Hi @Jannis Erz
Do you think it is feasible to work out something like this
for the approximation of
where
and then using the Symbolic Math Toolbox to transform it into a your desired Rational Polynomial function?
Jannis Erz
on 19 Sep 2022
Dear @Sam Chak
thanks a lot for the effort calculating the series expansion! However, I'm afraid that your proposed solution might be unfeasible for my application.
Let's discuss that down below!
Jannis Erz
on 19 Sep 2022
I tried to implement a first version regarding lsqcurvefit to a slightly different atan related function. However, I always get an error message as soon as I add rational parts to the approximate function. Could you be so kind and have a look to my code (see below)?
The rational function to approximate F_q has to look somehow like this one:
% Lateral Parameters
F_z0 = 3; % Nominal vertical wheel Load [kN]
F_z = 3000; % Actual vertical wheel load [kN]
muy0 = 1;
Cy = 1.3;
c1 = 8;
c2 = 1.33;
Ey = -1;
alpha_deg = [-20:0.5:20];
alpha = deg2rad(alpha_deg);
for idx1 = 1:length(alpha)
Dy0 = muy0 * (F_z/1000);
CFa0 = c1*c2*F_z0*sin(2*atan((F_z/1000)/(c2*F_z0))); % Cornering Stiffness [kN/rad]
By0 = CFa0/(Cy*Dy0);
df_z = ((F_z/1000)-F_z0)/F_z0; % Normalized change in vertical load % Normalized change in vertical load
F_q(idx1) = (Dy0 * sin(Cy*atan(By0*(alpha(idx1)) - Ey*(By0*(alpha(idx1))-atan(By0*(alpha(idx1)))))))*1000; % Pure Lateral Force [N]
end
fun = @(p,alpha)((p(1).*alpha.^3+p(2).*alpha.^2+p(3).*alpha)./((alpha.^4)+p(4).*alpha.^2)+p(5));
% Rational Function to approximate Lateral Tyre Force
p0 = [0.1,0.1,0.1,0.1,0.1];
p = lsqcurvefit(fun,p0,alpha,F_q);
plot(rad2deg(alpha),F_q,'o')
hold on
plot(rad2deg(alpha),fun(p,alpha))
xlabel('Tyre Side Slip Angle [deg]')
ylabel('Pure Lateral Tyre Force [N]')
legend('Pacejka Data','Fitted exponential')
Torsten
on 19 Sep 2022
fun = @(p,alpha)(p(1)*alpha.^3+p(2)*alpha.^2+p(3)*alpha)./(alpha.^4+p(4)*alpha.^2+p(5));
instead of
fun = @(p,alpha)((p(1).*alpha.^3+p(2).*alpha.^2+p(3).*alpha)./((alpha.^4)+p(4).*alpha.^2)+p(5));
Torsten
on 19 Sep 2022
Edited: Torsten
on 19 Sep 2022
The poles will make trouble.
Maybe you can formulate a constraint on p(4) and p(5) such that the roots are outside the interval of approximation. Then use fmincon instead of lsqcurvefit and define this constraint in the nonlcon function.
% Lateral Parameters
F_z0 = 3; % Nominal vertical wheel Load [kN]
F_z = 3000; % Actual vertical wheel load [kN]
muy0 = 1;
Cy = 1.3;
c1 = 8;
c2 = 1.33;
Ey = -1;
alpha_deg = [-20:0.5:20];
alpha = deg2rad(alpha_deg);
for idx1 = 1:length(alpha)
Dy0 = muy0 * (F_z/1000);
CFa0 = c1*c2*F_z0*sin(2*atan((F_z/1000)/(c2*F_z0))); % Cornering Stiffness [kN/rad]
By0 = CFa0/(Cy*Dy0);
df_z = ((F_z/1000)-F_z0)/F_z0; % Normalized change in vertical load % Normalized change in vertical load
F_q(idx1) = (Dy0 * sin(Cy*atan(By0*(alpha(idx1)) - Ey*(By0*(alpha(idx1))-atan(By0*(alpha(idx1)))))))*1000; % Pure Lateral Force [N]
end
fun = @(p)(p(1)*alpha.^3+p(2)*alpha.^2+p(3)*alpha)-F_q.*(alpha.^4+p(4)*alpha.^2+p(5));
fun1 = @(p)(p(1)*alpha.^3+p(2)*alpha.^2+p(3)*alpha)./(alpha.^4+p(4)*alpha.^2+p(5));
%fun = @(p,alpha)((p(1).*alpha.^3+p(2).*alpha.^2+p(3).*alpha)./((alpha.^4)+p(4).*alpha.^2)+p(5));
% Rational Function to approximate Lateral Tyre Force
p0 = [0.1,0.1,0.1,0.1,0.1];
format long
p = lsqnonlin(fun,p0)
Local minimum possible.
lsqnonlin stopped because the size of the current step is less than
the value of the step size tolerance.
p = 1×5
1.0e+03 *
1.573707413267325 0.000000000311668 -0.034491567121686 0.000043504411553 -0.000001348685846
rad2deg(roots([1 0 p(4) 0 p(5)]))
ans =
-0.000000000000003 +14.544020290876990i
-0.000000000000003 -14.544020290876990i
-8.289268225905060 + 0.000000000000000i
8.289268225905060 + 0.000000000000000i
plot(rad2deg(alpha),F_q,'o')
hold on
plot(rad2deg(alpha),fun1(p))
xlabel('Tyre Side Slip Angle [deg]')
ylabel('Pure Lateral Tyre Force [N]')
legend('Pacejka Data','Fitted exponential')
Torsten
on 19 Sep 2022
Edited: Torsten
on 19 Sep 2022
% Lateral Parameters
F_z0 = 3; % Nominal vertical wheel Load [kN]
F_z = 3000; % Actual vertical wheel load [kN]
muy0 = 1;
Cy = 1.3;
c1 = 8;
c2 = 1.33;
Ey = -1;
alpha_deg = [-20:0.1:20];
alpha = deg2rad(alpha_deg);
for idx1 = 1:length(alpha)
Dy0 = muy0 * (F_z/1000);
CFa0 = c1*c2*F_z0*sin(2*atan((F_z/1000)/(c2*F_z0))); % Cornering Stiffness [kN/rad]
By0 = CFa0/(Cy*Dy0);
df_z = ((F_z/1000)-F_z0)/F_z0; % Normalized change in vertical load % Normalized change in vertical load
F_q(idx1) = (Dy0 * sin(Cy*atan(By0*(alpha(idx1)) - Ey*(By0*(alpha(idx1))-atan(By0*(alpha(idx1)))))))*1000; % Pure Lateral Force [N]
end
fun = @(p)sum(((p(1)*alpha.^3+p(2)*alpha.^2+p(3)*alpha)./(alpha.^4+p(4)*alpha.^2+p(5))-F_q).^2);
fun1 = @(p)(p(1)*alpha.^3+p(2)*alpha.^2+p(3)*alpha)./(alpha.^4+p(4)*alpha.^2+p(5));
p0 = [0.1,0.1,0.1,0.1,0.1];
options = optimset('MaxFunEvals',100000,'MaxIter',100000);
format long
p = fmincon(fun,p0,[],[],[],[],[],[],@nonlcon,options)
Local minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.
p = 1×5
1.0e+03 *
1.320317528765901 0.000000004546610 0.059807026912929 0.000090641897411 0.000002053988392
rad2deg(roots([1 0 p(4) 0 p(5)]))
ans =
-0.000001716290581 +12.197536562689598i
-0.000001716290581 -12.197536562689598i
0.000001716290582 +12.197536562700607i
0.000001716290582 -12.197536562700607i
plot(rad2deg(alpha),F_q,'o')
hold on
plot(rad2deg(alpha),fun1(p))
xlabel('Tyre Side Slip Angle [deg]')
ylabel('Pure Lateral Tyre Force [N]')
legend('Pacejka Data','Fitted exponential')
function [c,ceq] = nonlcon(p)
c = p(4)^2/4-p(5);
ceq = [];
end
Jannis Erz
on 20 Sep 2022
Thank you so much for your effort handing me over such a good solution!
There is one thing still unclear to me regarding the nonlinear equality constraint of the denominator (x^4 + p4*x^2 + p5). I use x instead of alpha here for the sake of convenience. Why do we want the inner sqrt to become negative?
I do understand that the roots of the first solution are causing the serrated process output when the denominator gets close to zero, i.e., close to the root values (+- 8.3°) on the real axis.
Torsten
on 20 Sep 2022
One way to avoid real poles in the interval of interest is to make all roots of the denominator complex. That's what I did. Of course, this is quite arbitrary and might exclude a better solution. Feel free to experiment.
Jannis Erz
on 23 Sep 2022
@Torsten Again thanks a lot for the hints! I played a little with the roots and the best solutions seems to make all roots of the denominator complex. I also tried out the Curve Fitter App using the Rational Option (only available at R2022b). With one varying variable it delivers perfect fitting using rationals with the degree 3 and 4 respectively. Unfortunately it is limited to one variable only!
Now I'm trying to extend the problem to two varying variables (see code below). However, I'm still not satisfied with the results (absolute error around 500N). I'd like to force the absolute error below 200 N. @John D'Errico and @Torsten do you have any ideas how to systematically vary with the rational expressions and/or starting points to get better results?
%% Version 2 of Rational Approximation of tyre lateral force "F_q" with two varying variables side slip angle "a" and wheel load "F_z"
% Lateral Parameters
F_z0 = 3; % Nominal vertical wheel Load [kN]
muy0 = 1;
Cy = 1.3;
c1 = 8;
c2 = 1.33;
Ey = -1;
alpha_deg = -20:0.1:20;
a_vec = deg2rad(alpha_deg);
F_z_vec = linspace(500,7500,length(a_vec));
[a, F_z] = meshgrid(a_vec,F_z_vec);
for idx1 = 1:length(a)
for idx2 = 1:length(F_z)
Dy0 = muy0 * (F_z(idx1,idx2)/1000);
CFa0 = c1*c2*F_z0*sin(2*atan((F_z(idx1,idx2)/1000)/(c2*F_z0))); % Cornering Stiffness [kN/rad]
By0 = CFa0/(Cy*Dy0);
F_q(idx1,idx2) = (Dy0 * sin(Cy*atan(By0*(a(idx1,idx2)) - Ey*(By0*(a(idx1,idx2))-atan(By0*(a(idx1,idx2)))))))*1000; % Pure Lateral Force [N]
end
end
fun = @(p)sum((calcrational(p,a,F_z)-F_q).^2,'all'); % Problem to be minimized
fun1 = @(p)(calcrational(p,a,F_z)); % Rational polynomial approximation of F_q
p0(1,1:10) = 0.1;
options = optimset('MaxFunEvals',1000000,'MaxIter',1000000);
format long
p = fmincon(fun,p0,[],[],[],[],[],[],[],options)
F_q_approx = calcrational(p,a,F_z);
abs_err = abs(F_q_approx)-abs(F_q); % Absolute error
%mesh(a,F_z,F_q)
%hold on
%mesh(a,F_z,(calcrational(p,a,F_z)))
mesh(a,F_z,abs_err) % Display absolute error
xlabel('Side Slip Angle [rad]')
ylabel('Wheel Load [N]')
zlabel('Absolute Error [N]')
% Rational approximation function
function rational = calcrational(p,x,y)
num = p(1) + p(2)*x + p(3)*y + p(4)*x.^2 + p(5)*y.^2 + p(6)*x.*y;
denom = p(7) + p(8)*x + p(9)*y + p(10)*x.^2;
rational = num./denom;
end
Torsten
on 23 Sep 2022
Your problem is overfitted.
As you can see in MATLAB's rational models section, one of the parameters to be fitted must be normalized to 1.
Jannis Erz
on 23 Sep 2022
@Torsten what do you mean by normalized? Sorry, but I'm a little bit clueless here :)
How would I implement that in the code? Enforce a nonlinear equality constraint on, e.g., p(10) to be equal to 1?
Torsten
on 23 Sep 2022
Edited: Torsten
on 23 Sep 2022
No, in the same way as the rational polynomials in MATLAB are defined - by setting p(10) = 1 and fitting only 9 parameters.
Note that if p(1),...,p(10) is a solution for the Least-squares problem, then also p(1)*c,p(2)*c,...,p(10)*c is a solution for every constant c not equal to 0. This is prohibited by normalizing one of the parameters to 1.
Answers (1)
John D'Errico
on 18 Sep 2022
Edited: John D'Errico
on 18 Sep 2022
@Sam Chak has suggested an interesting idea in a comment. However, while it is not a bad idea at first glance, I suspect it will be problematic. The issue is, that classic atan series is not very strongly convergent. Expect to need many terms before you get anything good. And that means the polynomial approximation will be poor at best.
For example, how many terms do you need for convergence as a function of x to k significant digits, in the simple atan series?
For example, when x == 1, how many terms are required? We can look at it easily, since that is just the Leibniz formula for pi/4. Thus...
N = 0:1000;
S = mod(N+1,2)*2 - 1;
approx = cumsum(S./(2*N+1));
truth = pi/4;
semilogy(N,abs(approx - truth))
grid on
So we need thousands of terms for convergence near x==1.
However, if you can use range reduction methods to force the argument x to the atan to be in a small interval near zero, you get much faster convergence. The problem is, that in itself may take some work, and it is not at all trivial to get good convergence.
Instead, you may gain from using a direct rational polynomial approximation, perhaps something like those found in the classic, by JF Hart, et al, "Computer Approximations". (Start reading around page 120 in my edition from 1978. The tables there give some pretty good approxmations, though they still employ range reduction.)
Note: Even though Hart is an old text, it is still a book I love dearly, but that might apply only to a real gearhead numerical analyst like me. I think I recall it is available as a Dover reprint.
7 Comments
Sam Chak
on 19 Sep 2022
@John D'Errico, thanks for the heads-up. @Jannis Erz has mentioned about using the approximated Rational Polynomial function for some kind of Lyapunov Stability analysis.
Since the range for and are given, I guess that the stability analysis may be applied on the vicinity of some equilibrium points like and , and therefore
.
Does @Jannis Erz need to approximate until for the purpose of showing the stability proof around the equilibrium points?
Jannis Erz
on 19 Sep 2022
Dear @John D'Errico thanks a lot for the recommendations! I already got myself the 1978 version of your mentioned book! Let's see how I can get along with it :)
@Sam Chak You're right that I'm examining the equilibrium points for stability proof. In addition to that I'd like to analyse the domain of attraction. Therefore I need the approximation to be within -pi/4 and + pi/4.
John D'Errico
on 19 Sep 2022
Edited: John D'Errico
on 19 Sep 2022
@Sam Chak - the point is, there are arguably multiple ways to generate a rational polynomial approximation to a function. However, if you START with a poorly convergent Taylor series, you will likely not be happy with the result. This is my fear, and why I would suggest using a rational polynomial approximation itself for the atan which can have better properties over the desired domain.
@Jannis Erz - You can learn a lot from the Hart text. It remains one of the books I leave on my shelf next to my desk. (Abramowitz & Stegun is another, where my copy has tabs attached to the sections I have used regularly over the years.)
Bruno Luong
on 19 Sep 2022
Edited: Bruno Luong
on 19 Sep 2022
@John D'Errico: not necessary. It is well known that Padé fractional series converges better much than Taylor series, but the Padé series is derived directly from the Taylor series.
That being say, the best approximation on a given interval is just from numerical calculation.
Sam Chak
on 19 Sep 2022
OP has changed the function to:
F_z0 = 3;
F_z = 3000;
muy0 = 1;
Dy0 = muy0*(F_z/1000);
c1 = 8;
c2 = 1.33;
CFa0 = c1*c2*F_z0*sin(2*atan((F_z/1000)/(c2*F_z0)));
Cy = 1.3;
By0 = CFa0/(Cy*Dy0);
Ey = -1;
alpha_deg = linspace(-20, 20, 401);
alpha = deg2rad(alpha_deg);
Fq = 1000*Dy0*sin(Cy*atan(By0*alpha - Ey*(By0*alpha - atan(By0*alpha))));
plot(alpha_deg, Fq, 'linewidth', 1.5), ylim([-4e3 4e3])
grid on, xlabel('\alpha, deg'), ylabel('F_{q}')
Bruno Luong
on 19 Sep 2022
Yes, Padé fractional series is my loosly wording, rather use Padé approximant please.
Sam Chak
on 19 Sep 2022
help \sym\pade
PADE approximation of a symbolic expression
PADE(f <, x> <,options>) computes the third order Pade approximant of f at x = 0.
If x is not given, it is determined using symvar.
PADE(f, x, a <, options>) computes the third order Pade approximant of f at x = a.
A different order can be specified using the option 'Order' (see below).
The following options can be given as name-value pairs:
Parameter Value
'ExpansionPoint' a
Compute the Pade approximation about the point a. It is also
possible to specify the expansion point as third argument
without explicitly using a parameter value pair.
'Order' [m, n]
Compute the Pade approximation with numerator order m and
denominator order n.
'Order' m
Compute the Pade approximation with both numerator and
denominator order equal to m. By default, 3 is used.
'OrderMode' 'Absolute' or 'Relative'.
If the order mode is 'Relative' and f has a zero or pole at the
expansion point, add the multiplicity of that zero to the order.
If f has no zero or pole at the expansion point, this option has
no effect. Default is 'Absolute'.
Examples:
syms x
pade(sin(x))
returns (60*x - 7*x^3)/(3*(x^2 + 20))
pade(cos(x), x, 'Order', [1, 2])
returns 2/(x^2 + 2)
See also TAYLOR.
Documentation for pade
doc pade
See Also
Categories
Find more on Polynomials in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)