trying to find the polynomial equation representing the given four second-order differential equations
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=rhs(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=rhs(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=rhs(f2_Yb)
%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=rhs(f2_Ytheta)
%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 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.2ft x_f=%0.3ft y_f=%0.3ft x_r=%0.3ft y_r=%0.3ft x_b=%0.3ft y_b=%0.3ft x_theta=%0.3ft y_theta=%0.3fn’,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
%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=AB %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.5fn’,n,Cd)
%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
%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.^9For 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=rhs(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=rhs(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=rhs(f2_Yb)
%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=rhs(f2_Ytheta)
%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 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.2ft x_f=%0.3ft y_f=%0.3ft x_r=%0.3ft y_r=%0.3ft x_b=%0.3ft y_b=%0.3ft x_theta=%0.3ft y_theta=%0.3fn’,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
%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=AB %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.5fn’,n,Cd)
%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
%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 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=rhs(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=rhs(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=rhs(f2_Yb)
%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=rhs(f2_Ytheta)
%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 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.2ft x_f=%0.3ft y_f=%0.3ft x_r=%0.3ft y_r=%0.3ft x_b=%0.3ft y_b=%0.3ft x_theta=%0.3ft y_theta=%0.3fn’,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
%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=AB %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.5fn’,n,Cd)
%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
%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 solving multiple dependent second ord, numerical methods, polynomial equation representing ord MATLAB Answers — New Questions
