Azzera filtri
Azzera filtri

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

12 visualizzazioni (ultimi 30 giorni)
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
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)
%========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
end
  9 Commenti
Walter Roberson
Walter Roberson il 20 Lug 2024
% 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) = sym(x0);
yn_f = sym(y0);
xn_r(1) = sym(x0);
yn_r = sym(y0);
xn_b(1) = sym(x0);
yn_b = sym(y0);
xn_theta(1) = sym(x0);
yn_theta = sym(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=%s\t x_f=%s\t y_f=%s\t x_r=%s\t y_r=%s\t x_b=%s\t y_b=%s\t x_theta=%s\t y_theta=%s\n',string(tn(i)),string(xn_f(i)),string(yn_f(i)),string(xn_r(i)),string(yn_r(i)),string(xn_b(i)),string(yn_b(i)),string(xn_theta(i)),string(yn_theta(i)))
end
t=0 x_f=0 y_f=0 x_r=0 y_r=0 x_b=0 y_b=0 x_theta=0 y_theta=0 t=5 x_f=5*y_f y_f=(93750*Xb(t))/59 - (500*Y)/59 - (11593750*Xf(t))/59 - (271875*theta(t))/118 + (500*diff(Xb(t), t))/59 - (725*diff(theta(t), t))/59 x_r=5*y_r y_r=(12574*Xb(t))/9 - (100*Y)/9 - (28756*Xr(t))/9 + (873893*theta(t))/450 + (100*diff(Xb(t), t))/9 + (139*diff(theta(t), t))/9 x_b=5*y_b y_b=(5*exp(-5*t))/73 - (100*Y)/73 - (15662*Xb(t))/73 + (9375*Xf(t))/73 + (6287*Xr(t))/73 + (242741*theta(t))/3650 + (50*diff(Xf(t), t))/73 + (50*diff(Xr(t), t))/73 + (3*diff(theta(t), t))/73 x_theta=5*y_theta y_theta=(67*exp(-5*t))/2712 - (20173*Y)/13560 + (242741*Xb(t))/6780 - (90625*Xf(t))/904 + (873893*Xr(t))/13560 - (159290251*theta(t))/678000 + (5*diff(Xb(t), t))/226 - (725*diff(Xf(t), t))/1356 + (695*diff(Xr(t), t))/1356 t=10 x_f=10*y_f y_f=(187500*Xb(t))/59 - (1000*Y)/59 - (23187500*Xf(t))/59 - (271875*theta(t))/59 + (1000*diff(Xb(t), t))/59 - (1450*diff(theta(t), t))/59 x_r=10*y_r y_r=(25148*Xb(t))/9 - (200*Y)/9 - (57512*Xr(t))/9 + (873893*theta(t))/225 + (200*diff(Xb(t), t))/9 + (278*diff(theta(t), t))/9 x_b=10*y_b y_b=(10*exp(-5*t))/73 - (200*Y)/73 - (31324*Xb(t))/73 + (18750*Xf(t))/73 + (12574*Xr(t))/73 + (242741*theta(t))/1825 + (100*diff(Xf(t), t))/73 + (100*diff(Xr(t), t))/73 + (6*diff(theta(t), t))/73 x_theta=10*y_theta y_theta=(67*exp(-5*t))/1356 - (20173*Y)/6780 + (242741*Xb(t))/3390 - (90625*Xf(t))/452 + (873893*Xr(t))/6780 - (159290251*theta(t))/339000 + (5*diff(Xb(t), t))/113 - (725*diff(Xf(t), t))/678 + (695*diff(Xr(t), t))/678 t=15 x_f=15*y_f y_f=(281250*Xb(t))/59 - (1500*Y)/59 - (34781250*Xf(t))/59 - (815625*theta(t))/118 + (1500*diff(Xb(t), t))/59 - (2175*diff(theta(t), t))/59 x_r=15*y_r y_r=(12574*Xb(t))/3 - (100*Y)/3 - (28756*Xr(t))/3 + (873893*theta(t))/150 + (100*diff(Xb(t), t))/3 + (139*diff(theta(t), t))/3 x_b=15*y_b y_b=(15*exp(-5*t))/73 - (300*Y)/73 - (46986*Xb(t))/73 + (28125*Xf(t))/73 + (18861*Xr(t))/73 + (728223*theta(t))/3650 + (150*diff(Xf(t), t))/73 + (150*diff(Xr(t), t))/73 + (9*diff(theta(t), t))/73 x_theta=15*y_theta y_theta=(67*exp(-5*t))/904 - (20173*Y)/4520 + (242741*Xb(t))/2260 - (271875*Xf(t))/904 + (873893*Xr(t))/4520 - (159290251*theta(t))/226000 + (15*diff(Xb(t), t))/226 - (725*diff(Xf(t), t))/452 + (695*diff(Xr(t), t))/452 t=20 x_f=20*y_f y_f=(375000*Xb(t))/59 - (2000*Y)/59 - (46375000*Xf(t))/59 - (543750*theta(t))/59 + (2000*diff(Xb(t), t))/59 - (2900*diff(theta(t), t))/59 x_r=20*y_r y_r=(50296*Xb(t))/9 - (400*Y)/9 - (115024*Xr(t))/9 + (1747786*theta(t))/225 + (400*diff(Xb(t), t))/9 + (556*diff(theta(t), t))/9 x_b=20*y_b y_b=(20*exp(-5*t))/73 - (400*Y)/73 - (62648*Xb(t))/73 + (37500*Xf(t))/73 + (25148*Xr(t))/73 + (485482*theta(t))/1825 + (200*diff(Xf(t), t))/73 + (200*diff(Xr(t), t))/73 + (12*diff(theta(t), t))/73 x_theta=20*y_theta y_theta=(67*exp(-5*t))/678 - (20173*Y)/3390 + (242741*Xb(t))/1695 - (90625*Xf(t))/226 + (873893*Xr(t))/3390 - (159290251*theta(t))/169500 + (10*diff(Xb(t), t))/113 - (725*diff(Xf(t), t))/339 + (695*diff(Xr(t), t))/339
%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
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
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.426140e-22.
a1 = 5x1
27.3117 -0.7398 0.0371 -0.0003 0.0000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
%%%%%=========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)
27.3117 -0.7398 0.0371 -0.0003 0.0000
%========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)
Cd = 1.0000
fprintf('For n=%d. Coefficient of Determination=%0.5f\n',n,Cd)
For n=4. Coefficient of Determination=0.99998
%EQUATION FOR THE POLYNOMIAL
syms x y
%for n=2 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2
%for n=3 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2+a3*x.^3
%for n=4 activate lines below
a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);a4=a(:,5);
sympref("FloatingPointOutput",true);
y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4
y = 
%for n=5 activate lines below
%a0=a(:,1);a1=a(:,1);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
%for n=6 activate lines below
%a0=a(:,1);a1=a(:,1);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
%for n=7 activate lines below
%a0=a(:,1);a1=a(:,1);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
%for n=8 activate lines below
%a0=a(:,1);a1=a(:,1);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
%for n=9 activate lines below
%a0=a(:,1);a1=a(:,1);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

Accedi per commentare.

Risposte (0)

Prodotti


Release

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by