Я пытаюсь решить ode45

Question

Я пытаюсь решить проблему с ode45.

В ode45 используются x 'и x' '. то я хочу проехать от ode45 до уравнения x.

Я попытался решить с помощью функции polyfit.

polyfit(t,y,4) 

но произошла ошибка «векторы X и Y должны быть одного размера».

Я не знаю, что мне нужно.

По моему мнению, он может решить интегральную функцию. Это правильно?

Было бы очень признательно, если бы вы справились с этим и рассказали мне, как его решить.

код ниже:

мой код построен с основным сценарием и функции.

clear 
global L M m f Jw rw Fz Cd p A bw fw Nw g uw seta Te Tb y_curve X Q 
L=[0 0.025 0.05 0.1 0.125 0.15 0.175 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1]; 
M=[0 0.225 0.45 0.65 0.685 0.705 0.69 0.68 0.65 0.635 0.63 0.6275 0.625 0.6225 0.62 0.6175 0.615 0.6125 0.610 0.6075 0.6050 0.6 0.5975 0.5950]; 

m = 1400; f= 0.01; Jw = 0.65; rw= 0.31; Fz = 3560; Cd = 0.5; p = 1.202; 
A = 1.95; bw = 0; fw = 0; Nw = 4; g = 9.81; uw = 0; seta = 0; Te = 0; Tb = 1000; 
X=polyfit(L,M,20); 
x_curve = 0:0.025:1; 
y_curve = polyval(X,x_curve); 
[t,i] = ode45('dot',[0:0.1:1],[20 20]); 
Q = polyfit(t,i,4); 
%subplot(1, 2, 1); 
%plot(L,M,'p',x_curve,y_curve,'m'); 
xlabel('Lamda') 
ylabel('mu'); grid; 

%subplot(1, 2, 2); 
plot(t,i); grid on; 
xlabel('time(s)'); 
ylabel('x'); 

function xdot = dot(t,x) 
global m f Jw rw Fz Cd p A Nw Te Tb g Y X 
%UNTITLED2 Summary of this function goes here 
%Detailed explanation goes here 
xdot = zeros(2,1); 

lamda = (x(2)-x(1))/x(1); 

Y=X(1,1)*(lamda)^20+X(1,2)*(lamda)^19+X(1,3)*(lamda)^18+X(1,4)*(lamda)^17+X(1,5)*(lamda)^16+X(1,6)*(lamda)^15+X(1,7)*(lamda)^14+X(1,8)*(lamda)^13+X(1,9)*(lamda)^12+X(1,10)*(lamda)^11+X(1,11)*(lamda)^10+X(1,12)*(lamda)^9+X(1,13)*(lamda)^8+X(1,14)*(lamda)^7+X(1,15)*(lamda)^6+X(1,16)*(lamda)^5+X(1,17)*(lamda)^4+X(1,18)*(lamda)^3+X(1,19)*(lamda)^2+X(1,20)*(lamda)^1+X(1,21); 
xdot(1)=(-0.5*p*Cd*A*((x(1)*rw)^2)-f*m*g+(Nw*Y*Fz))/(rw*m); 
xdot(2)=(Te-Tb-rw*Y*Fz)/Jw; 
end 

Answer 1

Это сравнительно легко: выход i составляет 2 колонки в ширину (так как есть два элемента в x), но только один столбец в t. Таким образом:

Q1 = polyfit(t,i(:,1), 4); 
Q2 = polyfit(t,i(:,2), 4); 

делает трюк.

У меня, однако, также возникает проблема кондиционирования в вашей полиномиальной посадке. Эту проблему вы можете исправить, используя опцию scale/shift, как описано в help polyfit.

Существует также особенность где-то в вашей производной, так как последние несколько записей: NaN, вызванные одним из значений xdot без ограничений. Итак, вы получаете все NaN за Q1 и Q2.

Ниже приведен вариант кода, который реализует эти изменения, а также показать несколько приемов программирования, которые широко считается выше текущим стилем:

• Don't. Use. Global. Variables.
• Предпочитают функции над сценариями
• Документ о намерении каждой отдельной элементарной операции («блок команд»)
• Выбрать имена переменных, которые обеспечивают элементарную самодокументированию коды
• Не злоупотребляйте скобками. Используйте скобки, где необходимо, или когда они значительно улучшают чтение.
• Интервал между операторами математики. Сделать это намного проще, чтобы определить термины и понять, что принадлежит тому, что является фантастическим источником ошибок.
• Совместите свой код.Вы пишете код только один раз, но читайте его сотни раз - сделайте , что часть самая простая.
• В MATLAB: если вы выполняете повторяющееся задание на копирование или строку кода, превышающую отметку в столбце 85, возможно, это более короткий, более читаемый, более понятный и удобный для пользователя «векторизованный» способ сделать это.

Это работает здесь (R2010a):

function myFun 

    %// Your data 
    L = [0 0.025 0.05 0.1 0.125 0.15 0.175 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1]; 
    M = [0 0.225 0.45 0.65 0.685 0.705 0.69 0.68 0.65 0.635 0.63 0.6275 0.625 0.6225 0.62 0.6175 0.615 0.6125 0.610 0.6075 0.6050 0.6 0.5975 0.5950]; 

    %// Your coefficients 
    A = 1.95;  m = 1400; 
    bw = 0;  f = 0.01; 
    fw = 0;  Jw = 0.65; 
    Nw = 4;  rw = 0.31; 
    g = 9.81;  Fz = 3560; 
    uw = 0;  Cd = 0.5; 
    seta = 0;  p = 1.202; 
    Te = 0;  N = 20; %// (degree of polynomial) 
    Tb = 1000; 

    %// Find best-fitting polynomial 
    [X,~,mu] = polyfit(L,M,N); 

    %// Or, if you like, a much faster alternative: 
    %// X = bsxfun(@power, (L(:)-mean(L))/std(L), N:-1:0)\M(:); 

    %// Plot the poly over the data 
    %// DEBUG -- put a space between the percent sign and brace below to include 
    %{ 
    subplot(2,1,1), hold on 
    x_curve = 0:0.025:1; 
    y_curve = polyval(X,x_curve, [], mu); 
    plot(L,M,'p',x_curve,y_curve,'m'); 
    xlabel('Lambda') 
    ylabel('mu'); 
    grid on 
    %} 

    %// Carry out the integration 
    tspan = [0 1]; 
    x0 = [20 20]; 
    [t,i] = ode45(@dxdt, tspan, x0); 

    %// Fit a polynomal through the solutions 
    Q1 = polyfit(t,i(:,1), 4); 
    Q2 = polyfit(t,i(:,2), 4); 

    %// Plot integration results   
    %// DEBUG -- put a space between the percent sign and brace below to include 
    %{ 
    subplot(2,1,2), hold on 
    plot(t,i); 
    grid on  
    xlabel('time(s)'); 
    ylabel('x'); 
    %} 

    %// The time derivative 
    function xdot = dxdt(~,x) 

     %// Polynomial value for these X 
     Y = ((x(2)-x(1))/x(1)).^(N:-1:0) * X.'; %' 

     %// The derivative: 
     xdot = [ 
      (-p/2*Cd*A*(x(1)*rw)^2 - f*m*g + Nw*Y*Fz)/rw/m 
      (Te - Tb - rw*Y*Fz)/Jw 
     ]; 

    end  

end 

Обратите внимание, что %// и %' только там, чтобы сделать форматирование кода на SO работать для всех.