RK4 in MatLab Help!

Hey all,

I'm very new to MatLab. Our assignment is to solve an ODE using the RK4 Method. I'd appreciate some help on my code. Here's what I have so far. I am getting the error: Indexing cannot yield multiple results. Error in ==> RK4HW1 at 3 [t,v]=vexact(g,m,c_d,t);

Here is my code:
    %filename vexact.m
  %x = vexact( g,m,c_d,t,v, v0)
  %
  %Here is the vexact.m file:
  %function [t,v] = vexact(g,m,c_d,t)
  %
  %vexact = sqrt((g*m)/c_d)*tanh(sqrt((g*c_d)/m)*t);
  %
  %end

    [t,v]=vexact(g,m,c_d,t);
    t0=0;
    v=0;
    g=9.81;
    m=75;
    c_d=.25;
    v0=0;
    t_final = 15;
    dt = .3;
    N = ceil( (t_final - t0)/dt);
    t(1) = t0;
    v = v0;
    v(1) = v0;
    for i=1:N
        k1 = dt * feval(vexact,t,v);
        k2 = dt * feval(vexact, t+dt/2, v+k1/2);
        k3 = dt * feval(vexact, t+dt/2, v+k2/2);
        k4 = dt * feval(vexact, t+dt, v+k3);
        v = v + (k1+2*k2+2*k3+k4)/6;
        t = t + dt;
        T(i+1) = t;
        V(i+1,:) = v';
    end

I also need to plot the numerical & exact answers. I would greatly appreciate any help!

Answer 1

[Sean de Wolski]

So is this script also called vexact?

Also, when youc all feval, which is unnecessary and could be rewritten as

vexact(t,v)

your function vexact doesn't accept two of its inputs.

[Matt Tearle]

I was assuming that was just formatting problems in the question. But one way or another, it seems that MATLAB is seeing vexact as a variable, not a function, hence the indexing error.

Also: new dog?

[Sean de Wolski]

I realized the formatting after posting originally. Then discovered the recursion issues and undefined variables.

And, yes new dog, he retired from mushing this past winter and we got him at the beginning of the summer!

Answer 2

[Matt Tearle]

Ahhh the classic problem. A few things I notice:

1. Your vexact is, I assume, supposed to evaluate the exact solution. In that case, you shouldn't be evaluating vexact inside the RK4 solver. RK4 works by evaluating the rate equations (and numerically integrating them).
2. So you need to define a function that takes t and v (and parameters as needed) as inputs, and returns dv/dt.
3. Read up on how to define functions in MATLAB. Your definition line function [t,v] = vexact(g,m,c_d,t) says that the function name is vexact, the inputs are g, m, c_d, t and the outputs are t,v. This means you need to define variables t and v inside vexact.m. You define a variable vexact (bad idea, given that that's the name of the function). If t is an input (and you don't alter it), why have it as an output?
4. Preallocate your arrays T and V before the loop:

T = zeros(N+1,1);
V = zeros(N+1,length(v0));

[Matthew Alsheikh]

I'm brand new to this all, so I need all the help I can get!

Here's my new code:

function [t,v] = vexact(g,m,c_d,t)

g=9.81;

m=75;

c_d=.25;

t=0;

v=0;

vexact = sqrt((g*m)/c_d)*tanh(sqrt((g*c_d)/m)*t);

vstar = g - c_d/m * vexact * vexact;

v0=0;

t_final = 15;

dt = .3;

T=zeros(N+1,1);

V = zeros(N+1,length(v0));

N = ceil( (t_final - t)/dt);

t(1) = t;

v = v0;

v(1) = v0;

for i=1:N

k1 = dt * feval(vexact,t,v);

k2 = dt * feval(vexact, t+dt/2, v+k1/2);

k3 = dt * feval(vexact, t+dt/2, v+k2/2);

k4 = dt * feval(vexact, t+dt, v+k3);

v = v + (k1+2*k2+2*k3+k4)/6;

t = t + dt;

T(i+1) = t;

V(i+1,:) = v';

end

Any better?

[Matt Tearle]

What is the problem you're trying to solve? It looks suspiciously like a free-fall with air resistence. So, going by memory, the problem is: solve v' = g - c_m*v^2. The exact/analytic solution is v(t) = sqrt((g*m)/...).

So you need to:

1) write a function to evaluate v' as a function of t and v (with g and c_m as parameters)

2) write an RK4 solver

3) use your RK4 solver to get the numerical solution

4) evaluate the analytic solution

5) compare the answers from 3) & 4)

Note that you don't need a function file for the exact solution. You can make one if you want, but you don't have to. If you do, then you need two functions: one for vexact, one for vprime (or vdot or vstar or whatever you want to call it).

Also, the problem with your vexact function is that you're using vexact as the name of the function (function [t,v] = VEXACT(g,m,c_d,t)) and then later defining a variable with the same name (VEXACT = sqrt((g*m)/c_d)*tanh(sqrt((g*c_d)/m)*t);). Don't do this! The returns from your function (t and v) are defined by the declaration line, so you need to have t and v defined as variables by the time vexact.m is done executing.

While I'm here... don't use feval. Use a function handle. Once you've written a function that evaluates v' (= g - c_mv^2) called, for example, vprime.m, create a variable in your script that is a function handle to vprime: f = @vprime; Then to evaluate the function, you can just use f: k1 = dt*f(t,v);