obtaining single polynomial equation using Heun's method from 4 second order differential equations

Question

Jon 2024-8-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2142601-obtaining-single-polynomial-equation-using-heun-s-method-from-4-second-order-differential-equations

评论： Torsten 2024-8-4

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 个评论
显示 3更早的评论隐藏 3更早的评论

Jon 2024-8-4

For this project I am supposed to use 2 different numerical methods to find the polynomial. Although using "POLYFIT" would be easier that doen not fit what I am supposed to do. That is why I am tryong to find out how to solve the 4 dependent second-order differential equations. I thought uning "Function" but I get the errors above. It has also been suggested by a classmate to use "mymatlabFunction" but they didn't know how I could use it. That is what I am trying to figure out which would work best within the scope of the assignment.

Torsten 2024-8-4

在 MATLAB Online 中打开

    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.

请先登录，再进行评论。

请先登录，再回答此问题。