1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

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