obtaining single polynomial equation using Heun's method from 4 second order differential equations
6 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
I get this error when trying to solve these 4 second order differential equations using a program that solves them using Heun's method and then Gaussian elimination to get a polynomial equation the represents the system. Previsouly posted for help, but the soluitions have symbolics and I need a singular soluition that represents all 4 equations. Does any one have a suggestion on how I can slove this? Is using the "function" with all 4 equations the best way to do this, or is there a nicer method for doing this?
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]};
%-------------------------------------------------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
t0=0; %time at the start
tm=20; %what value of (t time) you are ending at
t=t0:tm; % time from time start to time finish 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
dx=5; %delta(x) or h
h=dx;
%==INPUT SECTION for Euler's and Heun's==%
fx=@(x,y,t)y;
fy=@(x,y,t)Projectfunction;
%==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=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=======%%%%%%%%
%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)
function [dydt] = Projectfunction(t,x,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
x1 = x(1); %x1=Xf
x2 = x(2); %x2=Xr
x3 = x(3); %x3=Xb
x4 = x(4); %x4=theta
y1 = y(1); %y1=Xf'=dXf/dt
y2 = y(2); %y2=Xr'=dXr/dt
y3 = y(3); %y3=Xb'=dXb/dt
y4 = y(4); %y4=theta'=dtheta/dt
dydt1 = (Ksf*((x3-(L1*x4)))-x1)+Bsf*((y3-(L1*y4)-y1)-(Kf*x1))/Mf;
dydt2 = (Ksr*((x3+(L2*x4))-x2)+Bsr*((y3+(L2*y4))-y2)-(Kr*x2))/Mr;
dydt3 = (Ksf*((x3-(L1*x4))-x1)+Bsf*((y3-(L1*y4))-y1)+Ksr*((x3+(L2*x4))-x2)+Bsr*((y3+(L2*y4))-y2)-(10*exp(-(5*t))))/Mb;
dydt4 = ((-(Ksf*(x1-(L1*x4))*L1)-(Bsf*(y1-(L1*y4))*L1)+(Ksr*(x2+(L2*x4))*L2)+(Bsr*(y2+(L2*y4))*L2)+((10*exp(-(5*t)))*L3)))/Ic;
[dydt] = [dydt1; dydt2; dydt3; dydt4];
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]};
5 Commenti
Torsten
il 4 Ago 2024
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;
This is not Heun's method.
It cannot happen that xn(i+1) and yn(i+1) appear on both sides of the updates because Heun's method is explicit.
Risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!