Error in assignment when solving RK2

i keep getting error saying error of assignment, how do i fix this?

clf;
W = 80; %Forcing wind

%initial condition
theta0 = 1e-3; %initial radians
z0 = [0;0;theta0;0];
tspan = [0, 1000]; %around 16 minutes

h = 1/1000; %using the frequency as step size
t = [tspan(1):h:tspan(2)];
n = length(t);

z = z0; Z = NaN(n,4); Z(1,:) = z';
   for i = 2:n
       z = rk2(@(t,y) fTacoma(t,y,W), t(i),h,z);
       Z(i,:) = z';
   end

     plot(t,Z);

function y = rk2(t,w,h)
k1 = zdot(t,z)
k2 = zdot(t+ (h/2) , w + (k1/2))

y = z + k2
return
end

Answer 1

[Walter Roberson]

Inside rk2 you have

k1 = zdot(t,z)

We do not know what zdot is; you have not defined any variable by that name, and there is no Mathworks function by that name.

If we presume that it is a function that you defined yourself, then:

you are trying to pass z to zdot. But inside rk2, z is not defined. You did define a variable named z outside of rk2, but you defined it in a script, and variables are never shared between scripts and functions that are called by the script.

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Luna]

yup zdot is another function

    function zdot = fTacoma(t,z, param) %ODE Tacoma Narrows
    W = param(1); %wind speed in km/h

    %Constants
    l = 6; %length across
    m = 2500; %mass in kg
    K = 1000; %spring constant in N
    d = 0.01; %damping coeeficient
    a = 0.2; %Hooke's non-linearlity coefficient
    w = 2*pi*(38/60); %frequency

    %ODE
    a1 = exp(a*(z(1)-l*sin(z(3))));
    a2 = exp(a*(z(1)+l*sin(z(3))));
    zdot(1,1) = z(2);
    zdot(2,1) = -d*z(2) - (K/(m*a))*(a1 + a2 -2) + 0.2*W*sin(w*t);
    zdot(3,1) = z(4);
    zdot(4,1) = -d*z(4) + (3*cos(z(3))/l)*(K/(m*a))*(a1 - a2);

    return
   end

am not sure how to fix it

[Walter Roberson]

That is not a function named zdot: that is a function named fTacoma.

Your line

         z = rk2(@(t,y) fTacoma(t,y,W), t(i),h,z);

invokes rk2 passing in four arguments:

1. the handle to the anonymous function @(t,y) fTacoma(t,y,W)
2. t(i)
3. h
4. z

but your function rk2 is

function y = rk2(t,w,h)

which expects t, w, h. The function handle would be passed in as t, t(i) would be passed in as w, h would be passed in as h, and the extra z passed in would trigger an error about too many parameters.

[Luna]

ok so i fixed up the code but not its just giving me endless answers

clf;
W = 80; %Forcing wind

%initial condition
theta0 = 1e-3; %initial radians
z0 = [0;0;theta0;0];
tspan = [0, 1000]; %around 16 minutes

h = 1/1000; %using the frequency as step size
t = [tspan(1):h:tspan(2)];
n = length(t);

z = z0; Z =  NaN(n,4); Z(1,:) = z';
   for i = 2:n
       z = rk2(@(t,z) fTacoma(t,z), t(i), z, h);
       Z(i,:) = z';
   end

     plot(t,Z)

function y = rk2(fTacoma, t, z, h)
k1 = h*fTacoma(t,z)
k2 = h*fTacoma(t+(h/2) , z+(k1/2))

y = z + k2
return
end

[Walter Roberson]

See enclosed