multiplying a function handle by a constant

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=======%%%%%%%%

[dpb]

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]

How do I do that? That is the question, and will it work if there are 4equations in the function handle?

[Torsten]

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]

As I said in the posting this works fine for a single 2nd ODE and for the project we are supposed to use this program as defined

[Torsten]

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]

ans = 8x1

1.0e+04 * 0.0005 0.0006 0.0007 0.0008 -4.0207 0.1141 0.0017 0.0248

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

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]

@Torsten so how would this be coded to work with 4-2nd ODE? The professor suggested using the matlabFunction for the differential equations, and never said that that line wouldn't work for multiple 2nd ODE equations. If he never said, maybe he didn't know, so then I need to understand how to fix this so that it would work.

[Torsten]

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]

@Torsten how would I use the output in Heun's? I could see using the "x" and "y", that you get from ode15s, into the least square method. The problem is the Professor wants us to use 2 numerical methods for the soluition and commands like ODE and Polyfit do not qualify so what I did is below how would the ode15s give me something I can use in Heun's?

 
%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
dx=5;                            %delta(x) or h
h=dx;


soln = ode15s(@fun,[t0 tm],[0 0 0 0 0 0 0 0])

%==INPUT SECTION for Euler's and Heun's==%
fx=@(x,y,t)y;
fy=@(x,y,t)matlabFunction(fun);

%==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 y=%0.3f\t z=%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=soln.x;
Y=soln.y;


%%__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=======%%%%%%%%



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);                    %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

% 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]};

Answer 1

[Steven Lord]

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)

Y = 1x5

1 4 9 16 25

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

twiceY = g(1:5)

twiceY = 1x5

2 8 18 32 50

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

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

ans = 1x5

2 2 2 2 2

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[dpb]

"...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.

[Steven Lord]

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)

y1 = 1x5

1 4 9 16 25

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

y2 = h(1:5)

y2 = 1x5

2 8 18 32 50

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

m = g(4);
y4 = m(1:5)

y4 = 1x5

4 16 36 64 100

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[dpb]

And... :)

Answer 2

[Star Strider]

‘... 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.

.

[Jon]

then how do I deal with in the for loop where it is adding and multiplying by a changing value? I get that I could create another function and multiply the original function by the constant and the new named function would deal with that, but wouldn't I run into the same problem then when adding? So for the following section of code

Eqns=[Xf_2; Xr_2; Xb_2; theta_2];
F1=matlabFunction(Eqns,'Vars',{[Xf Xr Xb theta],[Vf Vr Vb omega],'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

I get I can switch part of it to this

Eqns=[Xf_2; Xr_2; Xb_2; theta_2];
G1=matlabFunction(Eqns,'Vars',{[Xf Xr Xb theta],[Vf Vr Vb omega],'t'})
F1=@(x,y,t)h*F1(x,y,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));
    yn(i+1)=yn(i)+fy(xn(i),yn(i),tn(i));
    
   
    
    %==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

but wouldn't I run into the same issue when adding? so in the for loop how do I do this when the loop changes with each iteration?

for i=1:length(tn)
    %==EULER'S METHOD
    xn(i+1)=xn(i)+fx(xn(i),yn(i),tn(i));
    yn(i+1)=yn(i)+fy(xn(i),yn(i),tn(i));
    
    
    
    %==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

[dpb]

"...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]

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]

@dpb because if I understood you correctly, and I am not sure I did, I get the error "Conversion to double from function_handle is not possible"

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
H=[h h]



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];
G1=matlabFunction(Eqns,'Vars',{[Xf Xr Xb theta],[Vf Vr Vb omega],'t'})
F1=@(x,y,t)h*F1(x,y,t)
Exn=@(x,y,t,i)xn(i)+fx(xn(i),yn(i),tn(i))*h
Eyn=@(x,y,t,i)yn(i)+fy(xn(i),yn(i),tn(i))*h



%==INPUT SECTION for Euler's and Heun's==%
fx=@(x,y,t)y;
fy=@(x,y,t)G1;

%==CALCULATIONS SECTION==%
tn=t0:h:tm;
xn(1) = x0;
yn(1) = y0;


for i=1:length(tn)
    %==EULER'S METHOD
    xn(i+1)=Exn;
    yn(i+1)=Eyn;
    
   
    
    %==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=======%%%%%%%%

As you can see I added the Exn=@(x,y,t,i)xn(i)+fx(xn(i),yn(i),tn(i))*h

Eyn=@(x,y,t,i)yn(i)+fy(xn(i),yn(i),tn(i))*h

%==INPUT SECTION for Euler's and Heun's==%

fx=@(x,y,t)y;

fy=@(x,y,t)G1;

%==CALCULATIONS SECTION==%

tn=t0:h:tm;

xn(1) = x0;

yn(1) = y0;

for i=1:length(tn)

%==EULER'S METHOD

xn(i+1)=Exn;

yn(i+1)=Eyn;

but I dont think thats what you mean, as I said before I have never worked with function handles and I am learning this as I go so I'm struggling with understanding the behavior of the function handles.

[dpb]

dx=5;                              %delta(x) or h
h=dx;
...
F1=@(x,y,t)h*F1(x,y,t)

F1 = function_handle with value:
    @(x,y,t)h*F1(x,y,t)

F1(0,0,0)

Unrecognized function or variable 'F1'.

Error in solution>@(x,y,t)h*F1(x,y,t) (line 4)
F1=@(x,y,t)h*F1(x,y,t)

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]

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.

Answer 3

[Torsten]

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