1

我正在尝试解决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
4

1 回答 1

1

这相对容易：输出i是 2 列宽（因为 中有两个元素x），但t. 所以：

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

成功了。

但是，我在您的多项式拟合中也遇到了条件问题。您可以使用 scale/shift 选项解决此问题，如help polyfit.

您的导数中的某处也存在奇点，因为最后几个条目是由无限增长NaN的值之一引起的。xdot所以，你得到了所有NaN的Q1and Q2。

以下是实现此更改的代码版本，并展示了一些被广泛认为优于您当前风格的编程实践：

  • 不。利用。全球的。变量。
  • 更喜欢函数而不是脚本
  • 记录每个不同原子操作的意图（“命令块”）
  • 选择提供代码基本自我文档的变量名称
  • 不要过度使用括号。必要时使用括号，或者当它们显着提高阅读能力时使用。
  • 数学运算符之间的间距。更容易识别术语并理解什么属于什么——这是一个极好的错误来源。
  • 对齐你的代码。你只写一次代码，但读了几百次——让那部分变得最简单。
  • 在 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

请注意，%//and%'只是为了使 SO 上的代码格式适用于所有人。

于 2013-11-11T12:11:05.293 回答