data not showing up

IC=[0 0 0 0 0 0];
ti=0;
tf=20;
tspan=[ti,tf];
h = 1;

 %Scenario 1
fMu1=@(t,y)Mu(t,y,ti,tf,1);
[t1,y1]=ode45(fMu1,tspan,IC);
figure()
hold on
subplot(3,1,1)
plot(t1,y1(:,1),t1,y1(:,3),t1,y1(:,5))
legend('x1', 'x2', 'x3')
xlabel('t')
ylabel('positions');
title('Scenario 1: Positions')
subplot(3,1,2)
plot(t1,y1(:,2),t1,y1(:,4),t1,y1(:,6))
legend('v1', 'v2', 'v3')
xlabel('t')
ylabel('velocities');
title('scenario 1: velocities')
%excitation
%subplot(3,1,3)
hold off

 sc1res1 = myEuler3(fMu1,IC,ti,tf,h)
prob1res1 = array2table(sc1res1', 'VariableNames',{'t','x1','v1','x2','v2','x3','v3'})
sc1res2 = myHeun(fMu1,IC,ti,tf,h)
prob1res2 = array2table(sc1res2', 'VariableNames',{'t','x1','v1','x2','v2','x3','v3'})
sc1res3 = myRK4(fMu1,IC,ti,tf,h);
prob1res3 = array2table(sc1res3', 'VariableNames',{'t','x1','v1','x2','v2','x3','v3'})

When I run this function, the plots show up fine but in the myEuler3 and myHeun cases, all the columns except the first show up as all zeros. The myRK4 case also runs fine. I don't think it's a problem with the functions because they've run perfectly in other examples but i'll post them here anyway.

function [sol] = myEuler3(F, y0, ti, tf, h)
t = ti:h:tf;
y = zeros(length(y0),length(t));
y(:,1) = y0;
for i = 2:length(t)
y(:,i) = y(:,i-1) + h*F(t(i-1),y(:,i-1));
end
sol = [ t; y ];
end

function [sol] = myHeun(F, y0, ti, tf, h)
t = ti:h:tf;
y = zeros(length(y0),length(t));
y(:,1) = y0;
for i = 2:length(t)
  yp = y(:,i-1) + h*F(t(i-1),y(:,i-1));
  y(:,i) = y(:,i-1) + 0.5*h*(F(t(i-1),y(:,i-1))+F(t(i),yp));
end
sol = [ t; y ];
end

[Susan Santiago]

I forgot to add the Mu function

function f = Mu(t,y,ti,tf,n)
m1 = 55;
m2 = 400;
m3 = 100;
k1 = 230000;
k2 = 30000;
k3 = 50000;
k4 = 3000;
n2 = 1500;
n3 = 4000;
n4 = 700;
tm = (ti+tf)/2;
t1 = tf/4;
L0 = 5;
v = 15;
A = 0.03;

 M =[0 1 0 0 0 0;
    -(k1+k2)/m1 -(n2+n4)/m1 k2/m1 n4/m1 0 n4/m1;
    0 0 0 1 0 0;
    k2/m2 n2/m2 -(k3+k2)/m2 -(n3+n2)/m2 k3/m2 n3/m2;
    0 0 0 0 0 1;
    k4/m3 n4/m3 k3/m3 n3/m3 -(k3+k4)/m3 -(n3+n4)/m3];

 switch n
    case 1
    u = (A/2)*(1-cos(2*pi*v*t/L0));
    case 2
    u = 1*(t<tm)+(-1)*(t>=tm);
    case 3
    u = 0*(t<tm)+1*(t>=tm);
    case 4
    u = 1*(t<tm)+0*(t>=tm);
    case 5
    u = (A*t)*(t<=tm)+0*(t>tm);
    case 6
    u1 = A*t;
    u2 = A*(t-2*tm);
    u = u1.*(t<=tm)+u2.*(t>tm);
    case 7
    u1 =  A*t;
    u2 =  A*(-t+4*t1);
    u3 =  A*(-t+3*t1);
    u =   u1.*(t<=t1)+u2.*((t>=t1)&(t<=2*t1)) +u3.*(t>=2*t1);
    case 8
    u = A.*((t>=t1)&(t<=2*t1));
end

 vect=[0;-k1*u/m1;0;0;0;0];
f = M*y+vect;
end

[KALYAN ACHARJYA]

Undefined function or variable 'myRK4'? Can you clarify your question?

Which data not showing up, Is this third plot, %excitation?

[Susan Santiago]

as I said, the data for myEuler3 and myHeun are the ones with data not showing up. I didn't include the function for myRK4 since that one ran properly and I didn't want to put too many functions in the original problem but for clarity, here it is.

function [sol] = myRK4(dy, y0, ti, tf, h)
t = ti:h:tf;
y = zeros(length(y0),length(t));
y(:,1) = y0;
for i = 2:length(t)
t1 = t(i-1);
y1 = y(:,i-1);
k1 = dy(t1,y1);
t2 = t(i-1)+(0.5*h);
y2 = y(:,i-1)+(0.5*k1*h);
k2 = dy(t2,y2);
t3 = t(i-1)+0.5*h;
y3 = y(:,i-1)+0.5*k2*h;
k3 = dy(t3,y3);
t4 = t(i-1)+h;
y4 = y(:,i-1)+k3*h;
k4 = dy(t4,y4);
y(:,i) = y1 + (h/6)*(k1+2*k2+2*k3+k4); 
end
sol = [t ; y];
end

[Susan Santiago]

also the third plot, I haven't finished working on which is why it's commented

[KALYAN ACHARJYA]

Answer 1

[KALYAN ACHARJYA]

% Are you looking for this data

[Susan Santiago]

That's what shows up when I run it but that data is not correct. I don't know why i'm getting all zeros

[Stephan]

Hi Susan,

i couldnt solve your problem until now, but I was able to narrow the problem down a bit. The issue happens when you call the Input argument F of your function:

function [sol_euler3] = myEuler3(F, y0, ti, tf, h)
t = ti:h:tf;
y = zeros(length(y0),length(t));
y(:,1) = y0;
for i = 2:length(t)

% original:             y(:,i) = y(:,i-1) + h*F(t(i-1),y(:,i-1));
% works until here:     y(:,i) = y(:,i-1) +  h;
% doesnt work:          y(:,i) = F(t(i-1),y(:,i-1));

end
sol_euler3 = [ t; y ];
end

If you comment out the section "works until here" and execute the code, what you expect happens. However, invoking F (t (i-1), y (:, i-1)) always ends with the result 0, which always results in the result you see because of the multiplication and addition to the previous loop pass.

I'm not sure what you want to achieve when you call F. Anyway, here is your problem. That is also the reason why there is no problem in the function myRK4 - because F does not occur there.

Maybe you can describe what the call of F should actually reach / calculate (Should F (t, y) be the solution of your ode45 calculation?), then somebody here will be able to solve your problem quickly.

Best regards

Stephan

[Susan Santiago]

F is the function being evaluated. it's the same as the function being evaluated with ode45. When the function runs, the results from each method should be similar. How do I still solve using Euler's method while avoiding the issue that you point out?

[Stephan]

Hi Susan,

while trying to understand the issue i found that:

_

This is the result of your computation which is working. When i look at the dimensions you calculate i cant believe that this is a correct calculation.

The dimensions explode with each iteration step ending in a range from -2*10e116 up to 5*10e114 (!). Is these order of magnitude actually correct for your problem?

What im still trying to understand is, what F(t,y) is exactly doing (or should do...)

So far I have understood that you call the function Mu via the function handle fMu1, which you give t and y. Furthermore, I believe that you want to evaluate this calculation for a period of 20 seconds every second.

What I did not understand is where the y comes from, that you want to pass the function handle. In addition to the above-mentioned (suspected) problem, there will always be 0, because the function call at time 0 with the y0 = IC equals zero. In the next iteration step, you then evaluate from F (t = 1, y = 1) -> which corresponds to the call of the function of F (1, [0; 0; 0; 0; 0,0]) -> This yields again 0.

Best regards

Stephan