Solving Coupled System of ODEs in MATLAB

math_guy314

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?

swraman

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:

math_guy314

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.