Solving Coupled System of ODEs in MATLAB

  1. 1. The problem statement, all variables and given/known data

    I am asked to solve a coupled system of 5 ODEs. There is also a function, f, which describes the release of carbon dioxide over time. I am given the release rates at certain values of t and asked to interpolate for other values of t in the interval [1000 3000]. After solving, I am asked to plot f as a subplot on top of y(1).

    2. Relevant equations

    The differential equations are shown in the code below. In my MATLAB code, I use g instead of dy/dt. The other terms are contants, etc ... I don't really have a question about these. It's more about the MATLAB code itself.

    dy(1)/dt=(ps-y(1))/d+f/mu1;
    dy(2)/dt=((y(3)-y(2))*w-k1-(ps-y(1))*mu2/d)/vs;
    dy(3)/dt=(k1-(y(3)-y(2))*w)/vd;
    dy(4)/dt=((y(5)-y(4))*w-k2)/vs;
    dy(5)/dt=(k2-(y(5)-y(4))*w)/vd;

    3. The attempt at a solution

    function g = carbon(t,y)
    %constants
    d=8.64;
    mu1=495;
    mu2=.0495;
    vs=0.12;
    vd=1.23;
    w=.001;
    k1=.000219;
    k2=.0000612;
    k3=.997148;
    k4=.0679;
    %equations describing equilibrium between carbon dioxide and carbonate
    %dissolved in the shallow ocean
    h=(y(2)-((y(2))^2-k3*(y(4))*(2*(y(2))-y(4)))^.5)/k3;
    c=(y(4)-h)/2;
    ps=k4*h^2/c;
    %interpolation function for carbon dioxide rate
    yr=[1000 1850 1950 1980 2000 2050 2080 2100 2120 2150 2225 2300 2400 2500 3000];
    data=[0 0 1 4 5 8 10 10.5 10 8 3.5 2 0 0 0];
    f=interp1(yr,data,t,'pchip');
    %differential equations
    g=zeros(5,1);
    g(1)=(ps-y(1))/d+f/mu1;
    g(2)=((y(3)-y(2))*w-k1-(ps-y(1))*mu2/d)/vs;
    g(3)=(k1-(y(3)-y(2))*w)/vd;
    g(4)=((y(5)-y(4))*w-k2)/vs;
    g(5)=(k2-(y(5)-y(4))*w)/vd;
    end

    I tried typing this in the command window:
    %solution
    [T,Y]=ode15s(@carbon,[1000 3000],[1 2.01 2.23 2.2 2.26]);
    %graph
    subplot(2,1,1);
    plot(t,f)
    xlabel('Time, t', 'fontsize', 14)
    ylabel('Carbon Dioxide Release Rate, f(t)', 'fontsize', 14)
    subplot(2,1,2);
    plot(t,Y(:,1));
    xlabel('Time, t', 'fontsize', 14)
    ylabel('Partial Pressure of Carbon Dioxide, p(t)', 'fontsize', 14)

    The problem is that this is my first time using MATLAB and I am not getting my function to run properly. Can someone point me in the right direction?
     
  2. jcsd
  3. Assuming your differential equations are all typed in correctly, your ODE function looks ready to solve.

    Theres a couple things that I notice though:

    At the command window you run:

    >> plot(t,f)

    but there is no variable t or f in the workspace. Should the t be T?

    Also where is the f coming from? you never create f in your main workspace. Should this be the code you are trying to run:
    Code (Text):

    %solution
    [T,Y]=ode15s(@carbon,[1000 3000],[1 2.01 2.23 2.2 2.26]);
    %graph
    subplot(2,1,1);
    t = T
    yr=[1000 1850 1950 1980 2000 2050 2080 2100 2120 2150 2225 2300 2400 2500 3000];
    data=[0 0 1 4 5 8 10 10.5 10 8 3.5 2 0 0 0];
    f=interp1(yr,data,t,'pchip');

    plot(t,f)
    xlabel('Time, t', 'fontsize', 14)
    ylabel('Carbon Dioxide Release Rate, f(t)', 'fontsize', 14)
    subplot(2,1,2);
    plot(T,Y(:,1));
    xlabel('Time, t', 'fontsize', 14)
    ylabel('Partial Pressure of Carbon Dioxide, p(t)', 'fontsize', 14)
     
     
  4. Yes, thank you. That helps. I actually realized that I had t instead of T and fixed it. It took me a while to realize f didn't copy into the workspace. I think I understand now.
     
Know someone interested in this topic? Share this thead via email, Google+, Twitter, or Facebook

Have something to add?