Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: Runge kutta

  1. Aug 30, 2010 #1
    i would like to programing in matlab with runge kutta 4 with 2nd order the variables (x with y and Q) where x distance , y hight of water and Q discharge from the following tow equation

    dydx = (s0-(((Q0^2*n^2)/(B^2.*y^(10/3))*(1+2.*y/B)^(4/3))+dQdx))/(1-Q0^2/(g*B^2*y.^3));
    dQdx = -2*sqrt(2)/3*ce*(y-w).^1.5*Q0/(B^2*y.^2*sqrt(g));


    i traing the following program but some thing error can you correcting to me, with glad

    program 1

    function kutta4research2(x0,xe, y0,Q0,N)

    h = (xe-x0)/N; %the step size
    x(1) = x0;
    y(1) = y0; %the initial value
    Q(1)=Q0;
    for i = 1:N
    k1y = Q(i);
    k1Q = h*f(x(i), y(i), Q(i));
    k2y = Q(i)+k1Q*h/2;
    k2Q = h*f(x(i)+h/2, y(i)+(k1y)*h/2,Q(i)+k1Q*h/2);
    k3y = Q(i)+k2Q*h/2;
    k3Q = h*f(x(i)+h/2, y(i)+(k2y)*h/2,Q(i)+k2Q*h/2);
    k4y = Q(i)+k3Q*h;
    k4Q = h*f(x(i)+h, y(i)+k3y*h,Q(i)+k3Q*h);
    x(i+1)=x0+i*h;
    y(i+1) = y(i) + (k1y + 2*k2y + 2*k3y + k4y)*h/6;
    Q(i+1) = Q(i) + (k1Q + 2*k2Q + 2*k3Q + k4Q)*h/6;
    end

    [x' y' Q']

    %this line of code below is optional if visualization is not needed
    plot(x, y, 'r*')

    function dydx = f(x, y , Q)
    function dQdx = Q(i)

    y0=0.17;
    s0=0.0005;
    Q0=0.143;
    n=0.012;
    B=0.15;
    ce=0.785;
    w=0.15;
    g=9.81;
    dydx = (s0-(((Q0^2*n^2)/(B^2.*y^(10/3))*(1+2.*y/B)^(4/3))+dQdx))/(1-Q0^2/(g*B^2*y.^3));
    dQdx = -2*sqrt(2)/3*ce*(y-w).^1.5*Q0/(B^2*y.^2*sqrt(g));

    this program error


    program 2

    function kutta4research3(x0,xe, y0,Q0,N)

    h = (xe-x0)/N; %the step size
    x(1) = x0;
    y(1) = y0; %the initial value
    Q(1)=Q0;
    for i = 1:N
    k1y = Q(i);
    k1Q = h*f(x(i), y(i), Q(i));
    k2y = Q(i)+k1Q*h/2;
    k2Q = h*f(x(i)+h/2, y(i)+(k1y)*h/2,Q(i)+k1Q*h/2);
    k3y = Q(i)+k2Q*h/2;
    k3Q = h*f(x(i)+h/2, y(i)+(k2y)*h/2,Q(i)+k2Q*h/2);
    k4y = Q(i)+k3Q*h;
    k4Q = h*f(x(i)+h, y(i)+k3y*h,Q(i)+k3Q*h);
    x(i+1)=x0+i*h;
    y(i+1) = y(i) + (k1y + 2*k2y + 2*k3y + k4y)*h/6;
    Q(i+1) = Q(i) + (k1Q + 2*k2Q + 2*k3Q + k4Q)*h/6;
    end

    [x' y' Q']

    %this line of code below is optional if visualization is not needed
    plot(x, y,x,Q ,'r*')

    function dy = f(x, y , Q)

    y0=0.17;
    s0=0.0005;
    Q0=0.143;
    n=0.012;
    B=0.15;
    ce=0.785;
    w=0.15;
    g=9.81;
    th=12;
    dy = (s0-((((Q0^2*n^2)/B^2.*y^(10/3))*(1+2.*y/B)^(4/3))+2*sqrt(2)/3*ce*(y-w).^1.5*Q0/((B^2*y.^2*sqrt(g)))/(cos(th)-Q0^2/(g*B^2*y.^3))));


    this program Q value constant error



    program 3


    function [t,x,v] = rk4ode2(func, a, b, x0, xd0, h);
    % Solution of 2nd order ODE using Runge-Kutta 4th order
    % with constant step size. ODE solved is converted to
    % two 1st order equations. The RHS of the system is
    % dv/dt = func(t, x, v)
    % dx/dt = v
    % See for example rhs_smd.m for forced spring-mass-damper
    %
    % USAGE: [t, x, v] = rk4ode2(func,a,b,x0,xd0,h)
    %
    % input func = name of external function to evaluate the RHS
    % of the ODE (eg 'rhs_smd')
    % a, b = limits of integration
    % x0 = initial condition (position)
    % xd0 = initial condition (velocity)
    % h = stepsize
    %
    % output [t, x, v] = solution vectors

    t = [a];
    x = [x0];
    v = [xd0];
    i = 1;

    while t(i) < b

    k1x = v(i);
    k1v = feval(func, t(i) , x(i) , v(i) );

    k2x = v(i)+k1v*h/2;
    k2v = feval(func, t(i)+h/2, x(i)+k1x*h/2 , v(i)+k1v*h/2 );

    k3x = v(i)+k2v*h/2;
    k3v = feval(func, t(i)+h/2, x(i)+k2x*h/2 , v(i)+k2v*h/2 );

    k4x = v(i)+k3v*h;
    k4v = feval(func, t(i)+h , x(i)+k3x*h , v(i)+k3v*h );

    i = i+1;
    t(i) = t(i-1) + h;

    x(i) = x(i-1) + (k1x + 2*k2x + 2*k3x + k4x)*h/6;
    v(i) = v(i-1) + (k1v + 2*k2v + 2*k3v + k4v)*h/6;

    end

    this program i can not know how it run


    and if possible solution with finite difference and finite element
     
  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted