(adsbygoogle = window.adsbygoogle || []).push({}); 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?

**Physics Forums - The Fusion of Science and Community**

# Solving Coupled System of ODEs in MATLAB

Know someone interested in this topic? Share a link to this question via email,
Google+,
Twitter, or
Facebook

- Similar discussions for: Solving Coupled System of ODEs in MATLAB

Loading...

**Physics Forums - The Fusion of Science and Community**