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.
multiplying a function handle by a constant
5 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
Torsten
2024 年 8 月 6 日
編集済み: Torsten
2024 年 8 月 6 日
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).
回答 (3 件)
Steven Lord
2024 年 8 月 5 日
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 件のコメント
Steven Lord
2024 年 8 月 5 日
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)
Star Strider
2024 年 8 月 5 日
編集済み: Star Strider
2024 年 8 月 5 日
‘... 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 件のコメント
dpb
2024 年 8 月 5 日
編集済み: dpb
2024 年 8 月 6 日
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
2024 年 8 月 6 日
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.
Torsten
2024 年 8 月 6 日
編集済み: Torsten
2024 年 8 月 6 日
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
0 件のコメント
参考
カテゴリ
Help Center および File Exchange で Number Theory についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
また、以下のリストから Web サイトを選択することもできます。
南北アメリカ
- América Latina (Español)
- Canada (English)
- United States (English)
ヨーロッパ
- 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)
アジア太平洋地域
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)