Integrate inside ode45, every value of a a vector

Hi everyone,

I have the following part inside a script

if t<=1.5

x=(0:2.51326406981556e-05:2*pi)'; %limits for integration

.

.

.

for i=1:1:250001

Y=(((r.*cos(x)-y(5))-r.*cos(x)+y(6))./((((((r.*cos(x)-y(5))-r.*cos(x)+y(6))).^2+((r.*sin(x)-y(5))+r.*sin(x)).^2).^0.5))); % r is difined in the beginning of the script

Integral(i)=trapz(x,Y);

end

.

.

.

y(5)=y(6);

dy(6)=(-kr*(y(5)+L*y(11))+b*(-kr*(y(6)+L*y(12)))-0.36*Fload*Integral(i))/mcd; % Fload is defined in the beginning.

I want after every time step, which is 1e-5, to use the values of y(5) and y(6) in order to calculate the Integral(i) and use this value for the next calculation of y(5) and y(6).

Practically, i want to integrate every value of Y, which is a function of y(5), y(6), inside [0,2pi], and y(5) and y(6) to be calculated based on this integral.

Is this the correct way for doing this?

Thank you very much.

Answer 1

[darova]

no need of for loop

if t<=1.5
    x=(0:2.51326406981556e-05:2*pi)'; %limits for integration
    Y=(((r.*cos(x)-y(5))-r.*cos(x)+y(6))./((((((r.*cos(x)-y(5))-r.*cos(x)+y(6))).^2+((r.*sin(x)-y(5))+r.*sin(x)).^2).^0.5))); % r is difined in the beginning of the script
%     C1 = r*cos(x) - y(5) + y(6) - r*cos(x);
%     C2 = 2*r.*sin(x)-y(5);
%     Y=C1./sqrt(C1.^2+C2.^2); % r is difined in the beginning of the script
    IntegralY=trapz(x,Y);    
end

Show the equations you use

[Ilias Minas]

Unfortunately when i run it without for loop, it calculates only the first value of y vector and it gives NaN for the rest of the calculations.

I am solving the following set of differential equations with ode45. IntegralY is participating at dy(5) only.

Without the integral it gives correct results.

Thank you

if t<=1.5

u=t/3;

Fload = -16.16429+86.16429*exp(7.07658*u);

kcs = (609.7485*exp(7.07658*u)).*2000;

x=(0:2.51326406981556e-05:2*pi)'; %limits for integration

% Y=(((r.*cos(x)-y(5))-r.*cos(x)+y(6))./((((((r.*cos(x)-y(5))-r.*cos(x)+y(6))).^2+((r.*sin(x)-y(5))+r.*sin(x)).^2).^0.5)));

C1 = r*cos(x) - y(5) + y(6) - r*cos(x);

C2 = 2*r.*sin(x)-y(5);

Y=C1./sqrt(C1.^2+C2.^2);

IntegralY=trapz(x,Y);

dy(1)=y(2);

dy(2)=(-kcs*(y(1)-y(9)+r*y(3)-r*y(11)-r*y(7))+b*(-kcs*(y(2)-y(10)+r*y(4)-r*y(12)-r*y(8)))-Fload)/mpp ;

dy(3)=y(4);

dy(4) = (-kcs*(y(1)-y(9)+r*y(3)-r*y(11)-r*y(7))*r -ktiltpp*y(3)+b*(-kcs*(y(2)-y(10)+r*y(4)-r*y(12)-r*y(8))*r)-Fload*r)/Ippx;

dy(5)=y(6);

dy(6)=(-kr*(y(5)+L*y(11))+b*(-kr*(y(6)+L*y(12)))-0.36*Fload*IntegralY)/mcd;

dy(7)=y(8);

dy(8)= (kcs*r*y(1)+(r^2)*kcs*y(3)-y(7)*(2*r^2*kcs+ktiltcd)-2*kcs*r*y(9)-2*(r^2)*kcs*y(11)+b*(kcs*r*y(2)+(r^2)*kcs*y(4)-y(8)*(2*r^2*kcs+ktiltcd)-2*kcs*r*y(10)-2*(r^2)*kcs*y(12)))/Icd;

dy(9)=y(10);

dy(10)= (kcs*y(1)+r*kcs*y(3)-y(7)*(2*r*kcs)-2*kcs*y(9)-2*r*kcs*y(11)+b*(kcs*y(2)+r*kcs*y(4)-y(8)*(2*r*kcs)-2*kcs*y(10)-2*r*kcs*y(12)))/mcd;

dy(11)=y(12);

dy(12)= (r*kcs*y(1)+(r^2)*kcs*y(3)-kr*L*y(5)-y(7)*(2*(r^2)*kcs)-2*kcs*r*y(9)-(2*(r^2)*kcs+kin*L^2+kr*L^2)*y(11)+b*(r*kcs*y(2)+(r^2)*kcs*y(4)-kr*L*y(6)-y(8)*(2*(r^2)*kcs)-2*kcs*r*y(10)-(2*(r^2)*kcs+kin*L^2+kr*L^2)*y(12)))/Iinsh;

[darova]

Show the whole code

[Ilias Minas]

function dy = transient_dampednew(t,y)

dy = zeros(12,1);

mpp = 1.85;

Ippx = 0.014;

mcd = 1.41;

Icd = 0.0049;

Iinsh = 0.0064;

L = 0.109;

r = 0.125;

ktiltcd = 40000;

ktiltpp = 70000;

kin = 139988679;

kr = 7000000;

b = 2.5E-6;

f t<=1.5

u=t/3;

Fload = -16.16429+86.16429*exp(7.07658*u);

kcs = (609.7485*exp(7.07658*u)).*2000;

x=(0:2.51326406981556e-05:2*pi)'; %limits for integration

% Y=(((r.*cos(x)-y(5))-r.*cos(x)+y(6))./((((((r.*cos(x)-y(5))-r.*cos(x)+y(6))).^2+((r.*sin(x)-y(5))+r.*sin(x)).^2).^0.5)));

C1 = r*cos(x) - y(5) + y(6) - r*cos(x);

C2 = 2*r.*sin(x)-y(5);

Y=C1./sqrt(C1.^2+C2.^2);

IntegralY=trapz(x,Y);

dy(1)=y(2);

dy(2)=(-kcs*(y(1)-y(9)+r*y(3)-r*y(11)-r*y(7))+b*(-kcs*(y(2)-y(10)+r*y(4)-r*y(12)-r*y(8)))-Fload)/mpp ;

dy(3)=y(4);

dy(4) = (-kcs*(y(1)-y(9)+r*y(3)-r*y(11)-r*y(7))*r -ktiltpp*y(3)+b*(-kcs*(y(2)-y(10)+r*y(4)-r*y(12)-r*y(8))*r)-Fload*r)/Ippx;

dy(5)=y(6);

dy(6)=(-kr*(y(5)+L*y(11))+b*(-kr*(y(6)+L*y(12)))-0.36*Fload*IntegralY)/mcd;

dy(7)=y(8);

dy(8)= (kcs*r*y(1)+(r^2)*kcs*y(3)-y(7)*(2*r^2*kcs+ktiltcd)-2*kcs*r*y(9)-2*(r^2)*kcs*y(11)+b*(kcs*r*y(2)+(r^2)*kcs*y(4)-y(8)*(2*r^2*kcs+ktiltcd)-2*kcs*r*y(10)-2*(r^2)*kcs*y(12)))/Icd;

dy(9)=y(10);

dy(10)= (kcs*y(1)+r*kcs*y(3)-y(7)*(2*r*kcs)-2*kcs*y(9)-2*r*kcs*y(11)+b*(kcs*y(2)+r*kcs*y(4)-y(8)*(2*r*kcs)-2*kcs*y(10)-2*r*kcs*y(12)))/mcd;

dy(11)=y(12);

dy(12)= (r*kcs*y(1)+(r^2)*kcs*y(3)-kr*L*y(5)-y(7)*(2*(r^2)*kcs)-2*kcs*r*y(9)-(2*(r^2)*kcs+kin*L^2+kr*L^2)*y(11)+b*(r*kcs*y(2)+(r^2)*kcs*y(4)-kr*L*y(6)-y(8)*(2*(r^2)*kcs)-2*kcs*r*y(10)-(2*(r^2)*kcs+kin*L^2+kr*L^2)*y(12)))/Iinsh;

And it is solved using the following script

clear all

clc

t0=0;

tinterval=1e-5;

tend=1.5

zpp=0;

dzpp=0;

tpp=0;

dtpp=0;

ycd=0;

dycd=0;

tcd=0;

dtcd=0;

zcd=0;

dzcd=0;

tinsh=0;

dtinsh=0;

tspan =[t0:tinterval:tend];

initcond = [zpp dzpp tpp dtpp ycd dycd tcd dtcd zcd dzcd tinsh dtinsh];

options = odeset ('RelTol', 1e-5,'AbsTol',[1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5]);

[t,y]=ode45('transient_dampednew',tspan,initcond,options);

[darova]

Where is end of if statement?

[Ilias Minas]

There are two end below the last equation...

One for the function and one of the if statement.

Accidentaly i didnt paste it

[darova]

Can't figure out where those NaN come from. Just take number without NaN and calculate

ix = ~isnan(Y);
IntegralY=trapz(x(ix),Y(ix));

Can you show you original equations? Why do you need to integrate each step?

[Ilias Minas]

I have this equation in which the integral is the IntegralY function in the code. Υcd is the y(5) and Vcd is y(6).

I want the integral between 0-2π at every value of Ycd and Vcd. Thats why i need the integral at every time step.

[darova]

Everything looks correct

[Ilias Minas]

Thats what i am taking...strange

[darova]

Here is what i have for tend=1.5e-02

dashed line: using trapz

solid line: integral=1