Solving Coupled System of ODEs in MATLAB

In summary, the conversation revolved around solving a coupled system of 5 ODEs and creating a function to describe the release of carbon dioxide over time. The person was asked to solve and plot the function in MATLAB, but was having trouble with the code. After some troubleshooting, the code was fixed by correcting variable names and creating the missing variable f in the main workspace.
  • #1
math_guy314
2
0

Homework Statement



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).

Homework 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;

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 news on Phys.org
  • #2
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:
%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)
 
  • #3
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.
 

Related to Solving Coupled System of ODEs in MATLAB

1. How do I define a coupled system of ODEs in MATLAB?

To define a coupled system of ODEs in MATLAB, you need to first create a function that represents the derivatives of each variable in your system. This function should take in the current values of the variables as input and return the corresponding derivatives. Then, you can use the built-in ODE solving functions, such as ode45 or ode23, to solve the system.

2. What is the difference between ode45 and ode23?

The main difference between ode45 and ode23 is the order of accuracy. ode45 uses a variable step-size Runge-Kutta method, which is more accurate but may take longer to compute. ode23, on the other hand, uses a fixed step-size method, which is less accurate but faster.

3. How do I specify initial conditions for my system of ODEs?

You can specify initial conditions for your system by providing an array of initial values to the ODE solving function. The order of the values in the array should match the order of the variables in your ODE function. Alternatively, you can also use the setInitialConditions function to specify initial conditions after defining your ODE system.

4. Can I solve a system of ODEs with non-constant coefficients?

Yes, you can solve a system of ODEs with non-constant coefficients in MATLAB. You will need to define your ODE function to take in not only the current values of the variables, but also the current time as input. Then, you can use a for loop to solve the system at different time points, updating the coefficients as needed.

5. How can I plot the solutions of my coupled ODE system in MATLAB?

To plot the solutions of your system, you can use the plot function in MATLAB. You will need to first solve your system using one of the ODE solving functions, and then use the output from the solver to plot the values of each variable over time. You can also use the subplot function to plot multiple variables on the same graph for comparison.

Similar threads

  • Engineering and Comp Sci Homework Help
Replies
2
Views
972
  • Engineering and Comp Sci Homework Help
Replies
4
Views
1K
  • Engineering and Comp Sci Homework Help
Replies
1
Views
1K
  • Engineering and Comp Sci Homework Help
Replies
7
Views
1K
  • Engineering and Comp Sci Homework Help
Replies
1
Views
1K
  • Engineering and Comp Sci Homework Help
Replies
1
Views
1K
  • Engineering and Comp Sci Homework Help
Replies
1
Views
1K
  • Engineering and Comp Sci Homework Help
Replies
1
Views
2K
  • Engineering and Comp Sci Homework Help
Replies
3
Views
2K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
5
Views
2K
Back
Top