multiplying a function handle by a constant
Show older comments
For the code below I am trying to find the polynomial equation that represents the system. There are 4-2nd ODE equations I have made them into 4-1st ODE the first order differentials become "y1,y2,y3and y4" and the 2nd order differentials become a first order. I then put all 4 into a matlabFunction, how do I multiple the function handle by the constant"h" with out getting the error "Operator '*' is not supported for operands of type 'function_handle'."?
The rest of the code works if working with a single differential equation fx(x,y,t) with only "x,y,t" but in my equations I have "Xf,Xr,Xb,theta" and "Vf,Vr,Vb,omega" that I have chosen to represent as "x1,x2,x3,x4" and "y1,y2,y3,y4" respectively. The next question is will I run into problems here as well?
I am not that expirenced working with matlabFunction command to know how to get this to work. The project requires me to use 2 different numerical methods to find the polynomial equation that best fits so I cannot use "Polyfit" to get the polynomial that best fits. Any suggestions that can help me to get this to work would be appreciated.
clear,close,clc
%______________________________________________________________SOLUITION_2_Heun's_Method_(for second order differential equations)_&__Least-Square_Nethod____________________________________________________________%
%4 Equations representing the system working with
% MfXf"=Ksf([Xb-(L1*theta)]-Xf)+Bsf([Xb'-(L1*theta')-Xf')-(Kf*Xf);
% MrXr"=Ksr([Xb+(L2*theta)]-Xr)+Bsr([Xb'+(L2*theta')]-Xr')-(Kr*Xr) ;
% MbXb"=Ksf([Xb-(L1*theta)]-Xf)+Bsf([Xb'-(L1*theta')]-Xf')+Ksr([Xb+(L2*theta)]-Xr)+Bsr([Xb'+(L2*theta')]-Xr')+fa(t);
% Ic*theta"={-[Ksf(Xf-(L1*theta))*L1]-[Bsf(Xf'-(L1*theta'))*L1]+[Ksr(Xr+(L2*theta))*L2]+[Bsr(Xr'+(L2*theta'))*L2]+[fa(t)*L3]};
clc
clear
%-------------------------------------------------SYSTEM_PARAMETERS----------------------------------------------------------------------------------------------------------------------------------------------------------
Ic=1356; %kg-m^2
Mb=730; %kg
Mf=59; %kg
Mr=45; %kg
Kf=23000; %N/m
Ksf=18750; %N/m
Kr=16182; %N/m
Ksr=12574; %N/m
Bsf=100; %N*s/m
Bsr=100; %N*s/m
L1=1.45; %m
L2=1.39; %m
L3=0.67; %m
t=[0:20]; % time from 0 to 20 seconds
n=4; %order of the polynomial
%Initial Conditions-system at rest therefore x(0)=0 dXf/dt(0)=0 dXr/dt(0)=0 dXb/dt(0)=0 dtheta/dt(0)=0 ;
%time from 0 to 20 h=dx=5;
x0=0; %x at initial condition
y0=0; %y at initial condition
t0=0; %t at the start
dx=5; %delta(x) or h
h=dx;
tm=20; %what value of (x) you are ending at
Xf = sym('x1'); %x1=Xf
Xr = sym('x2'); %x2=Xr
Xb = sym('x3'); %x3=Xb
theta = sym('x4'); %x4=theta
Vf = sym('y1'); %y1=Xf' = Vf = dXf/dt = Xf_1
Vr = sym('y2'); %y2=Xr' = Vr = dXr/dt = Xr_1
Vb = sym('y3'); %y3=Xb' = Vb = dXb/dt = Xb_1
omega = sym('y4'); %y4=theta'= omega = dtheta/dt = theta_1
t = sym('t');
% MfXf"=Ksf([Xb-(L1*theta)]-Xf)+Bsf([Xb'-(L1*theta')-Xf')-(Kf*Xf);
% Vf'=(1/Mf)(Ksf([Xb-(L1*theta)]-Xf)+Bsf([Xb'-(L1*theta')-Xf')-(Kf*Xf));
% MrXr"=Ksr([Xb+(L2*theta)]-Xr)+Bsr([Xb'+(L2*theta')]-Xr')-(Kr*Xr) ;
% Vr'=(1/Mr)(Ksr([Xb+(L2*theta)]-Xr)+Bsr([Xb'+(L2*theta')]-Xr')-(Kr*Xr)) ;
% MbXb"=Ksf([Xb-(L1*theta)]-Xf)+Bsf([Xb'-(L1*theta')]-Xf')+Ksr([Xb+(L2*theta)]-Xr)+Bsr([Xb'+(L2*theta')]-Xr')+fa(t);
% Vb'=(1/Mb)(-Ksf([Xb-(L1*theta)]-Xf)-Bsf([Xb'-(L1*theta')]-Xf')-Ksr([Xb+(L2*theta)]-Xr)-Bsr([Xb'+(L2*theta')]-Xr')+fa(t));
% Ic*theta"={-[Ksf(Xf-(L1*theta))*L1]-[Bsf(Xf'-(L1*theta'))*L1]+[Ksr(Xr+(L2*theta))*L2]+[Bsr(Xr'+(L2*theta'))*L2]+[fa(t)*L3]};
% omega'=(1/Ic){[-Ksf*Xf + Ksf((L1)^2)*theta)] + [-Bsf*Vf*L1 + Bsf*Vf*L1 + Bsf((L1)^2)*omega)] + [Ksr*Xr*L2 + Ksr((L2)^2)*theta)] + [Bsr*Vr*L2 + Bsr((L2)^2)*omega)*L2] + [fa(t)*L3]};
Xf_2=(Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Vb-(L1*omega)-Vf)-(Kf*Xf)))/Mf;
Xr_2=(Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Vb+(L2*omega))-Vr)-(Kr*Xr))/Mr;
Xb_2=(Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Vb-(L1*theta'))-Vf)+Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Vb+(L2*omega))-Vr)+(10*exp(-(5*t))))/Mb;
theta_2=((-(Ksf*(Xf-(L1^2*theta))*L1)-(Bsf*(Vf-(L1*omega))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*(Vr+(L2*omega))*L2)+((10*exp(-(5*t)))*L3)))/Ic;
Eqns=[Xf_2; Xr_2; Xb_2; theta_2];
F1=matlabFunction(Eqns,'Vars',{'x1','x2','x3','x4','y1','y2','y3','y4','t'})
%==INPUT SECTION for Euler's and Heun's==%
fx=@(x,y,t)y;
fy=@(x,y,t)F1;
%==CALCULATIONS SECTION==%
tn=t0:h:tm;
xn(1) = x0;
yn(1) = y0;
for i=1:length(tn)
%==EULER'S METHOD
xn(i+1)=xn(i)+fx(xn(i),yn(i),tn(i))*h;
yn(i+1)=yn(i)+fy(xn(i),yn(i),tn(i))*h;
%==NEXT 3 LINES ARE FOR HEUN'S METHOD
tn(i+1)=tn(i)+h;
xn(i+1)=xn(i)+0.5*(fx(xn(i),yn(i),tn(i))+fx(xn(i+1),yn(i+1),tn(i+1)))*h;
yn(i+1)=yn(i)+0.5*(fy(xn(i),yn(i),tn(i))+fy(xn(i+1),yn(i+1),tn(i+1)))*h;
fprintf('t=%0.2f\t x=%0.3f\t y=%0.3f\n',tn(i),xn(i),yn(i))
end
%%%LEAST SQUARE METHOD-FINDS POLYNOMIAL FOR GIVEN DATA SET%%%%%
%INPUT SECTION for Least-Square
X=xn;
Y=yn;
%%__CALCULATIONS SECTION__%%
k=length(X); %NUMBER OF AVAILABLE DATA POINTS
m=n+1; %SIZE OF THE COEFFICENT MATRIX
A=zeros(m,m); %COEFFICENT MATRIX
for j=1:m
for i=1:m
A(j,i)=sum(X.^(i+j-2));
end
end
B=zeros(m,1); %FORCING FUNCTION VECTOR
for i=1:m;
B(i)=sum(Y.*X.^(i-1));
end
a1=A\B %COEFFICIENTS FOR THE POLYNOMINAL--> y=a0+a1*x+a2*x^2....an*x^n CAN BE REPLACED BY GAUSSIAN ELIMINATION
%%%%%=========GAUSSIAN ELIMINATION TO FIND "a"========%%%%%%
%%%INPUT SECTION
%CALCULATION SECTION
AB=[A B]; %Augumentent matrix
R=size(AB,1); %# OF ROWS IN AB
C=size(AB,2); %# OF COLUMNS IN AB
%%%%FOWARD ELIMINATION SECTION
for J=1:R-1
[M,I]=max(abs(AB(J:R,J))); %M=MAXIMUM VALUE, I=LOCATION OF THE MAXIMUM VALUE IN THE 1ST ROW
temp=AB(J,:);
AB(J,:)=AB(I+(J-1),:);
AB(I+(J-1),:)=temp;
for i=(J+1):R;
if AB(i,J)~=0;
AB(i,:)=AB(i,:)-(AB(i,J)/AB(J,J))*AB(J,:);
end
end
end
%%%%BACKWARDS SUBSTITUTION
a(R)=AB(R,C)/AB(R,R);
for i=R-1:-1:1
a(i)=(AB(i,C)-AB(i,i+1:R)*a(i+1:R)')/AB(i,i);
end
disp(a)
syms X
P=0;
for i=1:m;
TT=a(i)*X^(i-1); %T=INDIVIDUAL POLYNOMIAL TERMS
P=P+TT;
end
display(P)
%========END OF GAUSSIAN ELIMINATION=======%%%%%%%%
8 Comments
dpb
on 5 Aug 2024
HINT: You can't multiply the handle itself by a constant, no, but you can multiply the expression that uses the handle...which is what you're trying to evauate.
Jon
on 5 Aug 2024
Torsten
on 5 Aug 2024
As noted before, your update formula for xn and yn must be wrong because xn(i+1) and yn(i+1) appear on both sides of the equation.
xn(i+1)=xn(i)+0.5*(fx(xn(i),yn(i),tn(i))+fx(xn(i+1),yn(i+1),tn(i+1)))*h;
yn(i+1)=yn(i)+0.5*(fy(xn(i),yn(i),tn(i))+fy(xn(i+1),yn(i+1),tn(i+1)))*h;
Jon
on 5 Aug 2024
For a single ODE, it will work. But not for a system of four 2nd order ODEs.
When you use this equation
xn(i+1)=xn(i)+0.5*(fx(xn(i),yn(i),tn(i))+fx(xn(i+1),yn(i+1),tn(i+1)))*h;
you change the xn(i+1) from Euler's method.
But in the equation
yn(i+1)=yn(i)+0.5*(fy(xn(i),yn(i),tn(i))+fy(xn(i+1),yn(i+1),tn(i+1)))*h;
the xn(i+1) from Euler's method has to be used, not the updated xn(i+1) from Heun's method you computed one line before.
Further, xn and yn must both be matrices of size (4 x nt) because you solve a system of four 2nd order ODEs. Another error in the assignment.
Why not simply writing a function that supplies the time derivatives of your 8 ordinary differential equations ?
fun(0,[1 2 3 4 5 6 7 8]) % Order of the input vector: [Xf,Xr,Xb,theta,Vf,Vr,Vb,omega]
function dydt = fun(t,y)
Ic=1356; %kg-m^2
Mb=730; %kg
Mf=59; %kg
Mr=45; %kg
Kf=23000; %N/m
Ksf=18750; %N/m
Kr=16182; %N/m
Ksr=12574; %N/m
Bsf=100; %N*s/m
Bsr=100; %N*s/m
L1=1.45; %m
L2=1.39; %m
L3=0.67; %m
Xf = y(1); %x1=Xf
Xr = y(2); %x2=Xr
Xb = y(3); %x3=Xb
theta = y(4); %x4=theta
Vf = y(5); %y1=Xf' = Vf = dXf/dt = Xf_1
Vr = y(6); %y2=Xr' = Vr = dXr/dt = Xr_1
Vb = y(7); %y3=Xb' = Vb = dXb/dt = Xb_1
omega = y(8); %y4=theta'= omega = dtheta/dt = theta_1
dydt = zeros(8,1);
dydt(1) = Vf;
dydt(2) = Vr;
dydt(3) = Vb;
dydt(4) = omega;
dydt(5) = (Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Vb-(L1*omega)-Vf)-(Kf*Xf)))/Mf;
dydt(6) = (Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Vb+(L2*omega))-Vr)-(Kr*Xr))/Mr;
dydt(7) = (Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Vb-(L1*theta'))-Vf)+Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Vb+(L2*omega))-Vr)+(10*exp(-(5*t))))/Mb;
dydt(8) = ((-(Ksf*(Xf-(L1^2*theta))*L1)-(Bsf*(Vf-(L1*omega))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*(Vr+(L2*omega))*L2)+((10*exp(-(5*t)))*L3)))/Ic;
end
Jon
on 6 Aug 2024
Use the above function "fun" and see whether calling it from one of the MATLAB ode integrators (e.g. ODE15s) gives reasonable results.
If you succeed, you can use it directly for Heun's method.
Note that you must supply eight initial values at t = 0:
[Xf(0),Xr(0),Xb(0),theta(0),Vf(0),Vr(0),Vb(0),omega(0)] (in this order).
Jon
on 6 Aug 2024
Answers (3)
You cannot directly multiply a function handle by a number. But you can multiple the value you receive by evaluating a function handle by a number, and if you do that inside an anonymous function you have a function handle that is like multiplying the original by a number.
f = @(x) x.^2; % function handle 1
g = @(x) 2*f(x); % evaluate f, multiply by a constant
Y = f(1:5)
twiceY = g(1:5)
This works with both anonymous functions and "named" function handles.
f2 = @sin;
g2 = @(x) 2*f2(x);
Y2 = f(1:5);
twiceY2 = g(1:5);
twiceY2./Y2 % all 2's
3 Comments
"...you can multipl[y] ... inside an anonymous function [so] you have a function handle that is like multiplying the original by a number."
f = @(x) x.^2; % function handle 1
g = @(x) 2*f(x); % evaluate f, multiply by a constant
And, if the constant is to be determined rather than a fixed value, then you can also put it inside the argument list...
f = @(x) x.^2; % function handle 1
g = @(x,c) c*f(x); % evaluate f, multiply by a variably-valued constant
NOTA BENE: If you define
f = @(x) x.^2;
g = @(x) c*f(x);
instead, the value c will have to have been defined prior to the definition of g and will be fixed inside the function handle g regardless of whether c is changed or not.
And if you want to make a function handle that makes function handles, that's possible too.
f = @(x) x.^2;
g = @(c) @(x) c*f(x); % Note that g doesn't accept x but ...
h = g(2); % the function handle _created by calling_ g does accept x
y1 = f(1:5)
y2 = h(1:5)
m = g(4);
y4 = m(1:5)
dpb
on 5 Aug 2024
And... :)
Star Strider
on 5 Aug 2024
Edited: Star Strider
on 5 Aug 2024
0 votes
‘... how do I multiple the function handle by the constant"h" with out getting the error "Operator '*' is not supported for operands of type 'function_handle'."?’
Evaluate the function with the appropriate arguments, then do the multiplication (and any other arithmetic operations).
NOTE —
I do not understand ‘fx’ and ‘fy’.
The definition for ‘fy’ should probably be:
fy=@(x,y,t)F1(x1,x2,x3,x4,y1,y2,y3,y4,t);
In your matlabFunction call, if you want ‘x1’...‘x4’ and likewise for ‘y1’...‘y4’ to be vectors instead of discrete variables, enclose them in square brackets:
F1=matlabFunction(Eqns,'Vars',{['x1','x2','x3','x4'],['y1','y2','y3','y4'],'t'})
They will lose their specific identities and will be labelled as ‘in1’ and ‘in2’ row vectors instead. (I did not run your code with those changes.)
EDIT — Added NOTE.
.
6 Comments
Jon
on 5 Aug 2024
Edited: Walter Roberson
on 5 Aug 2024
"...then how do I deal with in the for loop where it is adding and multiplying by a changing value? "
What's to keep you from doing as I showed in passing the values into the function? There's nothing limiting what you can do inside the anonymous function other than it has to fit on a single line--but if that's an issue, there's no rule that one MUST use an anonymous function--an m-file works just as well.
Walter Roberson
on 5 Aug 2024
G1=matlabFunction(Eqns,'Vars',{[Xf Xr Xb theta],[Vf Vr Vb omega],'t'})
F1=@(x,y,t)h*F1(x,y,t)
In order for that code to work, F1(x,y,t) has to be defined before the second line is executed, in which case you are replacing F1 with a multiple of a call to the "shadowed" version of F1 that it had before. It would seem to make more sense if the call was
F1=@(x,y,t)h*G1(x,y,t)
Also,
fy=@(x,y,t)F1;
defines fy as a function handle that expects three inputs, and ignores those inputs, and returns the function handle F1 instead.
yn(i+1)=yn(i)+fy(xn(i),yn(i),tn(i));
That is going to fail, because fy returns a function handle, and you cannot add anything to a function handle.
Jon
on 5 Aug 2024
dx=5; %delta(x) or h
h=dx;
...
F1=@(x,y,t)h*F1(x,y,t)
F1(0,0,0)
You have defined F1 in terms of F1 -- building the anonymous function can't tell, but when you try to evaluate it, then it fails because F1 isn't an argument so it isn't available that way, nor was it already defined in the scope of the code that created the anonymous function so it can't be embedded into the function as implicit constant. Whatever F1 is supposed to that is being scaled by h, has to already be defined there or be an argument, one.
Exn=@(x,y,t,i)xn(i)+fx(xn(i),yn(i),tn(i))*h
Has an issue as discussed above; xn, yn, tn are NOT the same as the dummy arguments x, y, t -- the arrays xn, yn, tn at the time this code line is executed will be embedded as implicit values in the anonymous function and they will not change no matter what happens to them outside the function...just as h above and here as well. The function as written will be totally independent of whatever values you pass as x, y, t and will only return the ith value of the expression based on the embedded arrays.
If fx is
fx=@(x,y,t)y;
still at this point, there's no point in making it a function; it is equivalent to whatever y is passed and there's no point at all in having the other arguments since the expression isn't dependent upon them in any form. If y is an array, it will return the array, if it's a single value, it will return it.
It's not clear what your thought process was at the time you wrote that so not at all certain what it might have intended to be/do...
Eyn=@(x,y,t,i)yn(i)+fy(xn(i),yn(i),tn(i))*h
Walter Roberson
on 6 Aug 2024
G1=matlabFunction(Eqns,'Vars',{[Xf Xr Xb theta],[Vf Vr Vb omega],'t'})
fy=@(x,y,t)G1;
That defines fy as a function handle that accepts three parameters. When invoked, it ignores all three parameters and returns the function handle G1.
You need to invoke G1, such as
fy=@(x,y,t)G1(x,y,t)
But the short form would just be
fy=G1
which copies the function handle G1 into fy, leaving fy as a function handle that accepts three parameters and calculates the result of Eqns on the parameters.
This is Heun's method, but your solution explodes.
I get the same results with a MATLAB ODE integrator like ODE45. You will have to check your differential equations.
tstart = 0.0;
tend = 2.0;
Y0 = [0,0,0,0,0,0,0,0];
h = 0.001;
T = tstart:h:tend;
[T,Y] = heun_method(@fun,T,Y0);
figure(1)
plot(T,Y)
[T,Y] = ode45(@fun,[tstart tend],Y0);
figure(2)
plot(T,Y)
function [T,Y] = heun_method(f,T,Y0)
T = T(:);
Y0 = Y0(:);
N = numel(T);
n = numel(Y0);
Y = zeros(n,N);
Y(:,1) = Y0;
for i = 2:N
t1 = T(i-1);
t2 = T(i);
y = Y(:,i-1);
h = t2 - t1;
f1 = f(t1,y);
ys = y + h*f1;
f2 = f(t2,ys);
Y(:,i) = y + 0.5*h*(f1+f2);
end
Y = Y.';
end
function dydt = fun(t,y)
Ic=1356; %kg-m^2
Mb=730; %kg
Mf=59; %kg
Mr=45; %kg
Kf=23000; %N/m
Ksf=18750; %N/m
Kr=16182; %N/m
Ksr=12574; %N/m
Bsf=100; %N*s/m
Bsr=100; %N*s/m
L1=1.45; %m
L2=1.39; %m
L3=0.67; %m
Xf = y(1); %y1=Xf
Xr = y(2); %y2=Xr
Xb = y(3); %y3=Xb
theta = y(4); %y4=theta
Vf = y(5); %y5=Xf' = Vf = dXf/dt = Xf_1
Vr = y(6); %y6=Xr' = Vr = dXr/dt = Xr_1
Vb = y(7); %y7=Xb' = Vb = dXb/dt = Xb_1
omega = y(8); %y8=theta'= omega = dtheta/dt = theta_1
dydt = zeros(8,1);
dydt(1) = Vf;
dydt(2) = Vr;
dydt(3) = Vb;
dydt(4) = omega;
dydt(5) = (Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Vb-(L1*omega)-Vf)-(Kf*Xf)))/Mf;
dydt(6) = (Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Vb+(L2*omega))-Vr)-(Kr*Xr))/Mr;
dydt(7) = (Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Vb-(L1*theta'))-Vf)+Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Vb+(L2*omega))-Vr)+(10*exp(-(5*t))))/Mb;
dydt(8) = ((-(Ksf*(Xf-(L1^2*theta))*L1)-(Bsf*(Vf-(L1*omega))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*(Vr+(L2*omega))*L2)+((10*exp(-(5*t)))*L3)))/Ic;
end
Categories
Find more on Numerical Integration and Differential Equations 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!
