ode45 doest not return some results as zeros

Hi there,

I trying to solve a set of differential equation using ode45 and OOP coding style

The problem is:

I have a 4-dof model, it supposed to return 4 column matrix result

But the first four columns are zeros. I have no idea why i that.

please see the attached files

I very appreciate you help

Thanks very much.

clc; close; clear all;

%General input
init_data = Init();

%Road input
init_road = InitRoadSingleBump();
road=RoadSingleBump(init_road);

%solving dif equations
solveeq = Solver(init_data, road);

t = linspace(0,2,1e3);

solveeq.solve(t)

classdef Init
    properties
        %số liệu chưa chuẩn
        mu = 200;
        mb = 2000;
        mst= 200;
        md = 80;
        
        ku = 870000;
        kb = 80000;
        kst= 2000;
        kd = 1000;
        
        cb = 8000;
        cst= 800;
        cd = 500;
    end
end
classdef InitRoadSingleBump
    properties
        d1 = 0.1;
        d2 = 0.05;
        v = 1/3.6;
    end
end
classdef RoadSingleBump
    properties
        initRoadSB
    end
    
    methods
        
        % Constructor
        function obj = RoadSingleBump(initRoadSBInput)
            obj.initRoadSB = initRoadSBInput;
        end
        
        function y = road(obj, t)
            d1 = obj.initRoadSB.d1;
            d2 = obj.initRoadSB.d2;
            v = obj.initRoadSB.v;
            y = d2.*(sin(pi*v.*t/d1)).^2.*(t > 0 & t < d1/v) + 0*(t <= 0);
        end
        
        function plotting(obj,t)
           plot(t,obj.road(t));
        end
    
    end
end
classdef Solver    
    properties
        init_data;
        road;
    end
    
    methods
        function obj=Solver(init_data,road)
        obj.init_data=init_data;
        obj.road=road;
        end
        
        function solve(obj, t)
            [T,X]  = ode45(@obj.PTVP, t, [0; 0; 0; 0; 0; 0; 0; 0]);  
            plot(T, X);
        end
        
        function dxdt = PTVP(obj, t, x)
            y = obj.road.road(t);
            
            kb = obj.init_data.kb;
            mu = obj.init_data.mu;
            mb = obj.init_data.mb;
            mst= obj.init_data.mst;
            md = obj.init_data.md;
        
            ku = obj.init_data.ku;
            kst= obj.init_data.kst;
            kd = obj.init_data.kd;
        
            cb = obj.init_data.cb;
            cst= obj.init_data.cst;
            cd = obj.init_data.cd;
            
%     x(1) = vt mu
%     x(2) = vt mb
%     x(3) = vt mst
%     x(4) = vt md
%     x(5) = cv mu
%     x(6) = cv mb
%     x(7) = cv mst
%     x(8) = cv md

            dxdt(1) = x(1);
            dxdt(2) = x(2);
            dxdt(3) = x(3);
            dxdt(4) = x(4);
            dxdt(5) = (-ku *(x(5)-y)    +  kb *(x(6)-x(5)) +  cb *(x(2)-x(1))                  )/mu;
            dxdt(6) = (-kb *(x(6)-x(5)) +  kst*(x(7)-x(6)) -  cb *(x(2)-x(1)) + cst*(x(3)-x(2)))/mb;
            dxdt(7) = (-kst*(x(3)-x(2)) +  kd *(x(4)-x(3)) -  cst*(x(7)-x(6)) + cd *(x(8)-x(7)))/mst;
            dxdt(8) = (-kd *(x(4)-x(3)) -  cd *(x(8)-x(7))                                     )/md;
            dxdt    = dxdt(:);
            
        end
    end
end

Answer 1

[Jan]

The initial value is zero:

[T,X]  = ode45(@obj.PTVP, t, [0; 0; 0; 0; 0; 0; 0; 0]);  

The derivatives of the first 4 elements equal the value:

dxdt(1) = x(1);
dxdt(2) = x(2);
dxdt(3) = x(3);
dxdt(4) = x(4);

Then the value 0 in these components stays at 0, as expected.

[Torsten]

y is not zero ...

[Jan]

@Torsten: Correct. And therefore the components 5 to 8 differ from 0 also:

dxdt(5) = (-ku *(x(5)-y)    +  kb *(x(6)-x(5)) +  cb *(x(2)-x(1))                  )/mu;
dxdt(6) = (-kb *(x(6)-x(5)) +  kst*(x(7)-x(6)) -  cb *(x(2)-x(1)) + cst*(x(3)-x(2)))/mb;
dxdt(7) = (-kst*(x(3)-x(2)) +  kd *(x(4)-x(3)) -  cst*(x(7)-x(6)) + cd *(x(8)-x(7)))/mst;
dxdt(8) = (-kd *(x(4)-x(3)) -  cd *(x(8)-x(7)))/md;

This matches the explanation "But the first four columns are zeros", doesn't it?

[Torsten]

@Jan

Yes, you are right.

@Minh Tran

But I think the system itself is incorrect. Shouldn't it be

dxdt(1) = x(5);
dxdt(2) = x(6);
dxdt(3) = x(7);
dxdt(4) = x(8);

instead of

dxdt(1) = x(1);
dxdt(2) = x(2);
dxdt(3) = x(3);
dxdt(4) = x(4);

?

[Minh Tran]

ok, thanks very much guys, I'll correct it and inform you later on.

[Minh Tran]

it works now guys

            dxdt(1) = (-ku *(x(5)-y)    +  kb *(x(6)-x(5)) +  cb *(x(2)-x(1))                  )/mu;
            dxdt(2) = (-kb *(x(6)-x(5)) +  kst*(x(7)-x(6)) -  cb *(x(2)-x(1)) + cst*(x(3)-x(2)))/mb;
            dxdt(3) = (-kst*(x(3)-x(2)) +  kd *(x(4)-x(3)) -  cst*(x(7)-x(6)) + cd *(x(8)-x(7)))/mst;
            dxdt(4) = (-kd *(x(4)-x(3)) -  cd *(x(8)-x(7))                                     )/md;
            dxdt(5) = x(1);
            dxdt(6) = x(2);
            dxdt(7) = x(3);
            dxdt(8) = x(4);