Solving Coupled System of ODEs in MATLAB


by math_guy314
Tags: differential eqs, matlab, ode, plot
math_guy314
math_guy314 is offline
#1
Oct13-10, 01:04 PM
P: 2
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?
Phys.Org News Partner Science news on Phys.org
SensaBubble: It's a bubble, but not as we know it (w/ video)
The hemihelix: Scientists discover a new shape using rubber bands (w/ video)
Microbes provide insights into evolution of human language
swraman
swraman is offline
#2
Oct14-10, 07:50 AM
P: 145
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:
%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)
math_guy314
math_guy314 is offline
#3
Oct14-10, 09:08 PM
P: 2
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.


Register to reply

Related Discussions
Help needed for solving BVP for a system of ODEs Calculus & Beyond Homework 1
help in solving a system of ODEs using MATLAB Calculus & Beyond Homework 0
Solving ODEs system at Matlab Math & Science Software 17
Matlab- solving odes for a matrix Math & Science Software 0
Solving coupled ODEs with Matlab Math & Science Software 0