trying to find the polynomial equation representing the given four second-order differential equations

Question

0 votes

For the following problem in the first image i get the error pictured in the third image when trying to solve the 4 second-order differential equations using numerical methods. I am trying to use Heun's method for solving second-order differential equations and from here getting the x and y values corresponding to the function. Then taking these values and putting them into the least-square method to find the polynomial equation that best represents the system of equations. The third image is of the free-body diagram representing this system, then follows the code, if anyone can see where I need to adjust the code to run correctly this would be much appreciated.

MATLAB CODE:

% 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

%Initial Conditions x=(0)=0 dXf/dt(0)=y(0)=0 ;

%t from 0 to 1 h=dx=0.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

syms Xf(t) Xr(t) Xb(t) theta(t) fa(t) Y

Xf_1=diff(Xf,t);

Xf_2=diff(Xf,t,2);

Xr_1=diff(Xr,t);

Xr_2=diff(Xr,t,2);

Xb_1=diff(Xb,t);

Xb_2=diff(Xb,t,2);

theta_1=diff(theta,t);

theta_2=diff(theta,t,2);

%%%% MfXf"=Ksf((Xb-(L1*theta))-Xf)+Bsf((Xb'-(L1*theta'))-Xf')-(Kf*Xf);

%anything in [] below indicates a deritive

% Mf*[(d^2Xf)/(dt^2)] = Ksf*((Xb-(L1*theta))-Xf)+Bsf*(([(dXb)/(dt)]-(L1*[(dtheta)/(dt)])-[(dXf)/(dt)])-(Kf*Xf);

Mf*Xf_2-((Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Xb_1-(L1*theta_1)-Xf_1)-(Kf*Xf))))==0;

f2_YF=isolate(Mf*Xf_2-((Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Xb_1-(L1*theta_1)-Xf_1)-(Kf*Xf))))==0,Xf_2);

f2_Yf=subs(f2_YF,[Xf_2 Xf_1],[Xf_1 Y])

f2_Yf =

f2_Yf=rhs(f2_Yf)

f2_Yf =

%dXf/dt=y_f-->ASSUMED

%dXf/dt=f1_f(x,y,t)=y_f

%dYf/dt=f2_f(x,y,t)=f2_Yf

%%%% MrXr"=Ksr((Xb+(L2*theta))-Xr)+Bsr((Xb'+(L2*theta'))-Xr')-(Kr*Xr);

%anything in [] below indicates a deritive

% Mf*[(d^2Xr)/(dt^2)] = Ksf*((Xb-(L1*theta))-Xf)+Bsr*(([(dXb)/(dt)]-(L2*[(dtheta)/(dt)])-[(dXr)/(dt)])-(Kf*Xf);

Mr*Xr_2-((Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Xb_1+(L2*theta_1))-Xr_1)-(Kr*Xr)))==0;

f2_YR=isolate(Mr*Xr_2-((Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Xb_1+(L2*theta_1))-Xr_1)-(Kr*Xr)))==0,Xr_2);

f2_Yr=subs(f2_YR,[Xr_2 Xr_1],[Xr_1 Y])

f2_Yr =

f2_Yr=rhs(f2_Yr)

f2_Yr =

%dXr/dt=y_r-->ASSUMED

%dXr/dt=f1_r(x,y,t)=y_r

%dYr/dt=f2_r(x,y,t)=f2_Yr

%%%% MbXb"=Ksf([Xb-(L1*theta)]-Xf)+Bsf([Xb'-(L1*theta')]-Xf')+Ksr([Xb+(L2*theta)]-Xr)+Bsr([Xb'+(L2*theta')]-Xr')+fa(t);

%anything in [] below indicates a deritive fa(t) given as 10e^(-(5t))

% Mb*[(d^2Xb)/(dt^2)]=Ksf*((Xb-(L1*theta))-Xf)+Bsf*(([(dXb)/(dt)]-(L1*[(dtheta)/(dt)]))-[(dXf)/(dt)])+Ksr*([Xb+(L2*theta)]-Xr)+Bsr( (([(dXb)/(dt)]+(L2*[(dtheta)/(dt)]))-[(dXr)/(dt)]))+fa(t);

Mb*Xb_2-((Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Xb_1-(L1*theta_1))-Xf_1)+Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Xb_1+(L2*theta_1))-Xr_1))==(10*exp(-(5*t))));

f2_YB=Xb_2==(-Ksf*((Xb-(L1*theta))-Xf)-Bsf*((Xb_1-(L1*theta_1))-Xf_1)-Ksr*((Xb+(L2*theta))-Xr)-Bsr*((Xb_1+(L2*theta_1))-Xr_1)+(10*exp(-(5*t))))/Mb;

f2_Yb=subs(f2_YB,[Xb_2 Xb_1],[Xb_1 Y])

f2_Yb(t) =

f2_Yb=rhs(f2_Yb)

f2_Yb(t) =

%dXb/dt=y-->ASSUMED

%dXb/dt=f1_b(x,y,t)=y_b

%dYb/dt=f2_b(x,y,t)=f2_Yb

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

%anything in [] below indicates a deritive fa(t) given as 10e^(-(5t))

%Ic*[(d^2theta)/(dt^2)]=(-(Ksf*(Xf-(L1*theta))*L1)-(Bsf*([dXf/dt]-(L1*[(dtheta)/(dt)]))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*([(dXr)/(dt)]+(L2*[(dtheta)/(dt)]))*L2)+(fa(t)*L3));

theta_2==((-(Ksf*(Xf-(L1*theta))*L1)-(Bsf*(Xf_1-(L1*theta_1))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*(Xr_1+(L2*theta_1))*L2)+((10*exp(-(5*t)))*L3)))/Ic;

Ic*theta_2-(-(Ksf*(Xf-(L1*theta))*L1)-(Bsf*(Xf_1-(L1*theta_1))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*(Xr_1+(L2*theta_1))*L2))==(10*exp(-(5*t)))*L3;

f2_YTheta=theta_2==((-Ksf*(Xf-(Xb-(L1*theta)))*L1) - (Bsf*(Xf_1-(Xb_1-(L1*theta_1)))*L1) + (Ksr*(Xr-(Xb+(L2*theta)))*L2) + (Bsr*(Xr_1-(Xb_1+(L2*theta_1)))*L2) + ((10*exp(-(5*t)))*L3))/Ic;

f2_Ytheta=subs(f2_YTheta,[theta_2 theta_1],[theta_1 Y])

f2_Ytheta(t) =

f2_Ytheta=rhs(f2_Ytheta)

f2_Ytheta(t) =

%dXb/dt=y-->ASSUMED

%dXb/dt=f1_theta(x,y,t)=y_theta

%dYb/dt=f2_theta(x,y,t)=f2_Ytheta

syms x y_f y_r y_b y_theta t

%==INPUT SECTION==%

fx_f=@(x,y,t)y_f;

fy_f=@(x,y,t)f2_Yf;

fx_r=@(x,y,t)y_r;

fy_r=@(x,y,t)f2_Yr;

fx_b=@(x,y,t)y_b;

fy_b=@(x,y,t)f2_Yb;

fx_theta=@(x,y,t)y_theta;

fy_theta=@(x,y,t)f2_Ytheta;

%==CALCULATIONS SECTION==%

tn=t0:h:tm;

xn_f(1)=x0;

yn_f=y0;

xn_r(1)=x0;

yn_r=y0;

xn_b(1)=x0;

yn_b=y0;

xn_theta(1)=x0;

yn_theta=y0;

for i=1:length(tn)

%==EULER'S METHOD

xn_f(i+1)=xn_f(i)+fx_f(xn_f(i),yn_f(i),tn(i))*h;

yn_f(i+1)=yn_f(i)+fy_f(xn_f(i),yn_f(i),tn(i))*h;

xn_r(i+1)=xn_r(i)+fx_r(xn_r(i),yn_r(i),tn(i))*h;

yn_r(i+1)=yn_r(i)+fy_r(xn_r(i),yn_r(i),tn(i))*h;

xn_b(i+1)=xn_b(i)+fx_b(xn_b(i),yn_b(i),tn(i))*h;

yn_b(i+1)=yn_b(i)+fy_b(xn_b(i),yn_b(i),tn(i))*h;

xn_theta(i+1)=xn_theta(i)+fx_theta(xn_theta(i),yn_theta(i),tn(i))*h;

yn_theta(i+1)=yn_theta(i)+fy_theta(xn_theta(i),yn_theta(i),tn(i))*h;

%==NEXT 3 LINES ARE FOR HEUN'S METHOD

tn(i+1)=tn(i)+h;

xn_f(i+1)=xn_f(i)+0.5*(fx_f(xn_f(i),yn_f(i),tn(i))+fx_f(xn_f(i+1),yn_f(i+1),tn(i+1)))*h;

yn_f(i+1)=yn_f(i)+0.5*(fy_f(xn_f(i),yn_f(i),tn(i))+fy_f(xn_f(i+1),yn_f(i+1),tn(i+1)))*h;

xn_r(i+1)=xn_r(i)+0.5*(fx_r(xn_r(i),yn_r(i),tn(i))+fx_r(xn_r(i+1),yn_r(i+1),tn(i+1)))*h;

yn_r(i+1)=yn_r(i)+0.5*(fy_r(xn_r(i),yn_r(i),tn(i))+fy_r(xn_r(i+1),yn_r(i+1),tn(i+1)))*h;

xn_b(i+1)=xn_b(i)+0.5*(fx_b(xn_b(i),yn_b(i),tn(i))+fx_b(xn_b(i+1),yn_b(i+1),tn(i+1)))*h;

yn_b(i+1)=yn_b(i)+0.5*(fy_b(xn_b(i),yn_b(i),tn(i))+fy_b(xn_b(i+1),yn_b(i+1),tn(i+1)))*h;

xn_theta(i+1)=xn_theta(i)+0.5*(fx_theta(xn_theta(i),yn_theta(i),tn(i))+fx_theta(xn-theta(i+1),yn_theta(i+1),tn(i+1)))*h;

yn_theta(i+1)=yn_theta(i)+0.5*(fy_theta(xn_theta(i),yn_theta(i),tn(i))+fy_theta(xn-theta(i+1),yn_theta(i+1),tn(i+1)))*h;

fprintf('t=%0.2f\t x_f=%0.3f\t y_f=%0.3f\t x_r=%0.3f\t y_r=%0.3f\t x_b=%0.3f\t y_b=%0.3f\t x_theta=%0.3f\t y_theta=%0.3f\n',tn(i),xn_f(i),yn_f(i),xn_r(i),yn_r(i),xn_b(i),yn_b(i),xn_theta(i),yn_theta(i))

end

Unable to perform assignment because value of type 'sym' is not convertible to 'double'.

Caused by:
Error using mupadengine/feval2char
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.

%input

x=['x values from Heuns method'];

y=['y values from Heuns method'];

n=4; %order of the polynomial you want to find

%%%LEAST SQUARE METHOD-FINDS POLYNOMIAL FOR GIVEN DATA SET%%%%%

%INPUT SECTION

x=x; %x-cordinates from data-input; independent vairiables

y=y; %y-cordinates from data-output; dependent vairiables

%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

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

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

%========END OF GAUSSIAN ELIMINATION=======%%%%%%%%

%STANDARD DEVIATION

Y_bar=mean(y); %ADVERAGE OF y

St=sum((y-Y_bar).^2);

SD=sqrt(St/(k-1)); %STANDARD DEVIATION

%STANDARD ERROR

for i=1:m;

T(:,i)=a(i)*x.^(i-1); %T=INDIVIDUAL POLYNOMIAL TERMS

end

for i=1:k

y_hat(i)=sum(T(i,:));

end

Sr=sum((y-y_hat).^2);

Se=sqrt(Sr/(k-(n+1))); %STANDARD ERROR-Se

%COEFFICIENT OF DETERMINATION

Cd=(St-Sr)/St %COEFFICIENT OF DETERMINATION (r^2)

fprintf('For n=%d. Coefficient of Determination=%0.5f\n',n,Cd)

%EQUATION FOR THE POLYNOMIAL

syms x y

for i=0:n

for n=1

a0=a(:,1);a1=a(:,2);a2=a(:,3);

sympref("FloatingPointOutput",true);

y=a0+a1*x+a2*x.^2

end

for n=2

a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);

sympref("FloatingPointOutput",true);

y=a0+a1*x+a2*x.^2+a3*x.^3

end

for n=3

a0=a(:,1);a1=a(:,2);a2=a(:,2);a3=a(:,4);a4=a(:,5);

sympref("FloatingPointOutput",true);

y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4

end

for n=4

a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);

sympref("FloatingPointOutput",true);

y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5

end

for n=5

a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7)

sympref("FloatingPointOutput",true);

y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6

end

for n=6

a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);

sympref("FloatingPointOutput",true);

y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7

end

for n=7

a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);a8=a(:,9);

sympref("FloatingPointOutput",true);

y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7+a8*x.^8

end

for n=8

a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);a8=a(:,9);a9=a(:,10);

sympref("FloatingPointOutput",true);

y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7+a8*x.^8+a9*x.^9

end

for n=9

a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);a8=a(:,9);a9=a(:,10);a10=a(:,11);

sympref("FloatingPointOutput",true);

y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7+a8*x.^8+a9*x.^9+a10*x.^10

end

9 Comments
Show 7 older comments Hide 7 older comments

Jon on 19 Jul 2024

Edited: Torsten on 19 Jul 2024

Open in MATLAB Online

For the method that I am supposed to use he has us changing the first order into a symbolic y. I did forget to add the symbolic y_f y_r y_b and y_theta into the syms command line. the code after finding f2_Ytheta should appear as follows: