ODE45 for nonlinear problem

kindly check the below code....  

[madhan ravi]

Could you post the equation in LaTeX form?

[chiyaan chandan]

function code1_verification
t=0:0.05:50; % time scale
k1=100;
k2=50;
c1=0.2;
c2=0.4;
m1=10;
m2=20;
f1=@(t)(195*sin(t/2))/2 - 100*sin((4*t)/5)^3 + (5*cos(t/2) - 8*cos((4*t)/5))^3/5000;
f2=@(t)(8*cos((4*t)/5))/25 - 100*sin(t/2) - (16*sin((4*t)/5))/5 + 150*sin((4*t)/5)^3 - (5*cos(t/2) - 8*cos((4*t)/5))^3/5000;
[t,x]=ode45( @rhs, t, [0,0,0.5,0.8] );
 function dxdt=rhs(t,x)
 dxdt_1 = x(2);
 dxdt_2 = (f1(t)-c1*(x(2)-x(4)).^3-k1*(x(1)-x(3).^3))/m1;
 dxdt_3 = x(4);
 dxdt_4 = (f2(t)+c1*(x(2)-x(4)).^3 + k1*(x(1)-x(3).^3)+ k2*x(3).^3 - c2*x(4))/m2;
 dxdt=[dxdt_1; dxdt_2;dxdt_3;dxdt_4];
 end
clf
subplot(211),plot(t,x(:,1));
legend('displacement')
hold on 
subplot(212),plot(t,x(:,2));
legend('velocity')
xlabel('t'); ylabel('x');
end

--------------------------------------------------------------------------------------------------------------------------

i am not getting displacement and velocity profiles as sin(w*t) and w*cos(w*t) form respectively.

procedure

Two equations of motions are

m1*a1+c1*(v1-v2)^3+k1*(x1-x2^3)=f1..........(1)

m2*a2-c1*(v1-v2)^3-k1*(x1-x2^3)+k2*x2^3+c2*v2=f2..........(2)

i have considered

k1=100; k2=50; c1=0.2; c2=0.4; m1=10; m2=20;

displacement x1=sin(0.5*t)

velocity v1=0.5*cos(0.5*t)

Acceleration a1 =-0.5^2*sin(0.5*t);

displacement x2=sin(0.8*t)

velocity v2=0.8*cos(0.8*t)

Acceleration a2 = - 0.8^2*sin(0.8*t);

subtituted x1,v1,a1,x2,v2,a2 in above equation (1) and (2)

f1= (195*sin(t/2))/2 - 100*sin((4*t)/5)^3 + (5*cos(t/2) - 8*cos((4*t)/5))^3/5000

f2=8*cos((4*t)/5))/25 - 100*sin(t/2) - (16*sin((4*t)/5))/5 + 150*sin((4*t)/5)^3 - (5*cos(t/2)

Answer 1

[darova]

• avoid global. YOu can make nested function in your case

function linear_vs_nonlinear
% variables
% code
    function SD_L=linear(t,s)        
        % code
    end
end

• you have 4 equations, but only 2 initial conditions

• some operators are forgotten (multiplier i think)

[chiyaan chandan]

function code1_verification
t=0:0.05:50; % time scale
k1=100;
k2=50;
c1=0.2;
c2=0.4;
m1=10;
m2=20;
f1=@(t)(195*sin(t/2))/2 - 100*sin((4*t)/5)^3 + (5*cos(t/2) - 8*cos((4*t)/5))^3/5000;
f2=@(t)(8*cos((4*t)/5))/25 - 100*sin(t/2) - (16*sin((4*t)/5))/5 + 150*sin((4*t)/5)^3 - (5*cos(t/2) - 8*cos((4*t)/5))^3/5000;
[t,x]=ode45( @rhs, t, [0,0,0.5,0.8] );

 function dxdt=rhs(t,x)
 dxdt_1 = x(2);
 dxdt_2 = (f1(t)-c1*(x(2)-x(4)).^3-k1*(x(1)-x(3).^3))/m1;
 dxdt_3 = x(4);
 dxdt_4 = (f2(t)+c1*(x(2)-x(4)).^3 + k1*(x(1)-x(3).^3)+ k2*x(3).^3 - c2*x(4))/m2;
 dxdt=[dxdt_1; dxdt_2;dxdt_3;dxdt_4];
 end
clf
subplot(211),plot(t,x(:,1));
legend('displacement')
hold on 
subplot(212),plot(t,x(:,2));
legend('velocity')
xlabel('t'); ylabel('x');
end

--------------------------------------------------------------------------------------------------------------------------

i am not getting displacement and velocity profiles as sin(w*t) and w*cos(w*t) form respectively.

procedure

Two equations of motions are

m1*a1+c1*(v1-v2)^3+k1*(x1-x2^3)=f1..........(1)

m2*a2-c1*(v1-v2)^3-k1*(x1-x2^3)+k2*x2^3+c2*v2=f2..........(2)

i have considered

k1=100; k2=50; c1=0.2; c2=0.4; m1=10; m2=20;

displacement x1=sin(0.5*t)

velocity v1=0.5*cos(0.5*t)

Acceleration a1 =-0.5^2*sin(0.5*t);

displacement x2=sin(0.8*t)

velocity v2=0.8*cos(0.8*t)

Acceleration a2 = - 0.8^2*sin(0.8*t);

subtituted x1,v1,a1,x2,v2,a2 in above equation (1) and (2)

f1= (195*sin(t/2))/2 - 100*sin((4*t)/5)^3 + (5*cos(t/2) - 8*cos((4*t)/5))^3/5000

f2=8*cos((4*t)/5))/25 - 100*sin(t/2) - (16*sin((4*t)/5))/5 + 150*sin((4*t)/5)^3 - (5*cos(t/2)

[darova]

The code looks correct. I think you formulas are strange. Why displacement and velocity are in power 3?

I didn't understand what if f? Where is expression for it and why did you substituted x,v and a?

[chiyaan chandan]

f1 and f2 are the forcing terms.

i calculated f1 and f2 by substituting know displacement, velocity and acceleration terms. so that i can verify my dispalcement plot esily with sin(w*t), velocity with w*cos(w*t) and acceleration with -w^2*sin(w*t).

[chiyaan chandan]

kindly see the equations once.

[darova]

How do you know that ? Usually and are uknowns if you have ODE. and should be given

[chiyaan chandan]

To verify the results i assumed x1=sin(w1t), x2=sin(w2*t). Once the displcement, velocity and accelerations are obtained from the calculation. then that plot should show the sin(w1*t) form for diplacement, w1*cos(w1*t) for velocity. and -w1^2*sin(w1*t) for acceleration similarly for second system same form will be repeated by taking w2 value.

for this problem f1 and f2 are forcing terms are the time series input (acceleration vs time ) .before giving the that inputs. i need to check the results by appying the know forces. once it gives that as a output then i can go for the next step.