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 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
Torsten
il 6 Ago 2024
Modificato: Torsten
il 6 Ago 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).
Risposte (3)
Steven Lord
il 5 Ago 2024
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 Commenti
Steven Lord
il 5 Ago 2024
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
il 5 Ago 2024
Modificato: Star Strider
il 5 Ago 2024
‘... 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 Commenti
dpb
il 5 Ago 2024
Modificato: dpb
il 6 Ago 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
il 6 Ago 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.
Torsten
il 6 Ago 2024
Modificato: Torsten
il 6 Ago 2024
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 Commenti
Vedere anche
Categorie
Scopri di più su Number Theory in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
Puoi anche selezionare un sito web dal seguente elenco:
Americhe
- América Latina (Español)
- Canada (English)
- United States (English)
Europa
- 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-Pacifico
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)