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

Coupled first order differential equations

  1. Jul 1, 2005 #1
    How i can solve a system of 6 first order differential equations by using numerical techniques like Euler method, RK-4th order method , ODE -45 etc.How i can solve a system of 6 first order differential equations by using numerical techniques like Euler method, RK-4th order method , ODE -45 etc.
     
  2. jcsd
  3. Jul 3, 2005 #2
    can you give an example
     
  4. Jul 3, 2005 #3
    Coupled ODEs

    Dear here are the system of six differential equtions i want to solve them by Euler method , Predictor corrector method and for RK-4th order method. I have made a program in Matlab plz check wether it is right or not.
    function [GG ] = myfun_rafiq(t,x)

    % ----- State Variables Selection
    % x1 = nw, x2 = Nw, x3 = nc, x4 = Nc, x5 = np, x6 = Np

    % ----- Parameters

    sigmaEjQj = 793.7;
    EiQi =700;
    EpQp = 13.7;
    EcQc = 80;
    EfQf = 0;
    Vc = 9.08e6;
    Vw = 1.37e7;
    Vp = 1.37e6;
    lk = 0.0;
    Kc = 40.0;
    Kp = 6.90;
    f_t = 1;
    sigma = 13.4e-24;
    fn = 1;
    fs = 0.5;
    phi0 = 9.2e13;
    phiE = 0.0026*phi0;
    lambda = 7.4612e-005; % PER SECOND

    N0 = 6.023e23;
    A = 56;
    S = 1.01e8;

    % ---- parameters + definitio of C(t)

    a = 50*3600; % in seconds
    deltaCDeltaT = 20e-12;
    % deltaCDeltaT=deltaCDeltaT/(3600^2);
    C0 = 2.4e-13;
    Cs = 25e-6;
    b = 400*3600; % in seconds
    t0 =50*3600; % in seconds %REPLCED 500 with 50
    % t = 900*3600;
    if t < a
    C_t = 0;
    elseif (t > a & t < b)
    C_t = deltaCDeltaT*(t-t0);
    else
    C_t = Cs;
    end
    % t
    % C_t

    % ---- Parameters + definition of g(t)
    %
    tin = 500*3600;
    tmax = 600*3600;
    w0 = 18.3e6;
    w2 = 0.3*w0;
    alpha = 0.004;
    % g_t = 1;

    if t <= tin
    g_t = 1;
    elseif (t > tin & t <= tmax)
    g_t = 1 - alpha*(tin - t);
    else %REPLACED 'b' with tmax
    g_t = w2/w0;
    end
    % if(t<tin)
    % g_t = 1;
    % elseif((t<=tmax))
    % g_t = 1 - alpha*(tin - t);
    % % g_t=0.00001;
    % else %REPLACED 'b' with tmax
    % g_t = w2/w0;
    % % g_t=0.00001;
    % end

    % % pause

    Sw = (C_t* S*N0*fn*fs)/(Vw*A);

    dxdt = zeros(6,1);

    dxdt(1) = sigma*phiE*x(2) - ( (sigmaEjQj*g_t/Vw) + lk*g_t/Vw + lambda)*x(1) + (Kp*g_t/Vw)*x(5) + (Kc*g_t/Vw)*x(3);

    dxdt(2) = -( (sigmaEjQj*g_t/Vw) + lk*g_t/Vw + sigma*phiE)*x(2) + (Kp*g_t/Vw)*x(6) + (Kc*g_t/Vw)*x(4) + Sw;

    dxdt(3) = sigma*phi0*x(4) + (EcQc*g_t/Vc)*x(1) - ( Kc*g_t/Vc + lambda)*x(3);

    dxdt(4) = (EcQc*g_t/Vc)*x(2) - ( Kc*g_t/Vc + sigma*phi0)*x(4);

    dxdt(5) = (EpQp*g_t/Vp)*x(1) - ( Kp*g_t/Vp + lambda)*x(5);

    dxdt(6) = (EpQp*g_t/Vp)*x(2) - ( Kp*g_t/Vp)*x(6);

    GG=[dxdt(1); dxdt(2); dxdt(3); dxdt(4); dxdt(5); dxdt(6)];

    % ---- PIEAS

    % Clearing work space

    clc;
    clear all;
    % close all;
    tic
    % Define time for simulation
    %x0 = [0.1, 0.2, 0.3, 0.4, 0.1, 0.1]*10^-7;
    % x0=[0 0.5 0 0.5 0 0.5];
    % T = 1:500*3600; % 10 seconds
    % T=T*3600;

    % Define initional conditions

    %R=myfun_rahat(T,x0);

    % Run simulation
    % options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);

    [tt,xx] = ode45(@myfun_rafiq,[0 1000*3600],[0 0.5 0 0.5 0 50]);

    toc

    % Plot Results
    grid on;
    figure(1),plot( tt, xx(:,1),'r.-',tt,xx(:,3),'m:',tt,xx(:,5),'b.-');
    xlabel('t(sec)');
    title('Specific Activity');
    legend('nw','nc','np');
    figure(2)
    grid on;
    plot(tt, xx(:,1),'r-');
    %axis([0 1.81e5 0 3e-4]);
     
  5. Jul 4, 2005 #4

    saltydog

    User Avatar
    Science Advisor
    Homework Helper

    Oh my God dude! Now, I really admire you for taking the time to input all that code but it's just too hard to read. I picked out the first one above and I can't understand it. I realize you're probably not familiar with the math-formatting code we use in here called "LaTex" but that would help with understanding your code. If you want to try using it, you can go to the General Physics forum and read the "Introducing LaTex" thread at the top. Also, I use Mathematica so would not be able to help you with Matlab.
     
  6. Jul 4, 2005 #5

    HallsofIvy

    User Avatar
    Staff Emeritus
    Science Advisor

    If you are already familiar with numerical methods for a single equation in one variable, just repeat them: For 6 coupled equations in 6 unknowns, set up 6 Runge-Kutta iterations, doing all 6 at each step then using the results from all 6 for the values of the 6 variables in the next iteration.
     
  7. Jul 14, 2005 #6
    I had exactly the same problem when trying to simulate a double inverted pendulum a few months ago. The system is 3 coupled second order differential equations which I reduced to 6 first order differential equations. Here's the program (written in Delphi).
    http://atlas.walagata.com/w/peterbone/Balance.zip
    The code for the double inverted pendulum is in the DoublePendulum.pas file which you can open in notepad. You can then see how the equations are solved. You should be able to understand it even if you don't know pascal.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: Coupled first order differential equations
Loading...