MATLAB: How to graph a “system” of ODEs

differential equationsiterationode

I am trying to model a shape that is defined with multiple differential equations. For reference, the variables are x and y (from the normal Cartesian plane) and then S (arc length) and theta (angle of the arc length with the horizontal).

The equations are:

dtheta/dS = -sin(theta)/x + y – 2H (a constant)

dx/dS = cos(theta)

dy/dS = sin(theta)

I'm just wondering if MatLab has an easy way to do this. If not, how could I set up a modified Euler's Method to solve this?

Best Answer

Is this what you want to do?

          % dThXYdS(1) = theta(S), dThXYdS(2) = x(S),  dThXYdS(3) = y(S) 
          H = 10;
          dThXYdS = @(t,XYTh) [-sin(XThXY(1))/ThXY(2) + ThXY(3) - 2*H;  cos(ThXY(1));  sin(XYTh(1))];
          x0 = [0.1; 0.1; 0];
          Tspan = [0:0.01:2]';
          [T ThXY] = ode45(dThXYdS, Tspan, x0);
          figure(8)
          plot3(ThXY(:,1), ThXY(:,2), ThXY(:,3))
          xlabel('X(S)')
          ylabel('Y(S)')
          zlabel('\Theta(S)')
          grid

I have no idea what you intend for H or the rest, but this seems to work. I named the variable ThXY to denote the vector as [Theta(S); X(S); Y(S)].

NOTE: This is a duplicate of 'System of arc-length defined ODEs with ode45'.

Related Solutions

MATLAB: System of arc-length defined ODEs with ode45

As with your 'How can I graph a "system" of ODEs?' post, I suggest:

        % dThXYdS(1) = Theta, dThXYdS(2) = R(t) = x,  dThXYdS(3) = Z(t) 
        H = 10;
        dThXYdS = @(t,ThXY) [-sin(ThXY(1))/ThXY(2) + ThXY(3) - 2*H;  cos(ThXY(1));  sin(ThXY(1))];
        x0 = [0.1; 0.1; 0];
        Tspan = [0:0.01:2]';
        [T ThXY] = ode45(dThXYdS, Tspan, x0);
        figure(8)
        plot(ThXY(:,2), ThXY(:,3))
        xlabel('R(S)')
        ylabel('Z(S)')
        grid

Except for the axis labels in the plot, I used the variable designation from your previous post (and my previous answer) rather than change it to match your current variable designation.

The reason you are getting a matrix of NaNs is that your initial conditions are [0 0 0]. So R is zero when theta is zero and of course sin(theta) will be zero as well. By convention, (0/0) = NaN. An initial NaN in a recursive calculation such as yours creates a matrix of NaNs.

MATLAB: How to plot these 3 ode45 functions in 1 graph

The lines are not easily distinguishable.

Try this:

X0=2;
Sin=1000;
So=0;
S0=So+Sin;
Ks=150;
Y=0.5;
XS0=[2, 1000];
THv = [10 20 40];
for k = 1:numel(THv)
    tspan = [0 THv(k)];
    [t{k},XS{k}]=ode45(@(t,XS)BSfunction(t,XS,THv(k)), tspan, XS0);
end
figure
LineStl = {'-','-.','--'};
hold on
for k = 1:numel(THv)
    semilogy(t{k}, XS{k}, 'LineWidth',1.5, 'DisplayName', ['TH = ',num2str(THv(k))], 'LineStyle',LineStl{k})
end
hold off
grid
% lgdstr = compose('TH = %d', THv);
legend
figure
for k = 1:numel(THv)
    subplot(numel(THv), 1, k)
    plot(t{k}, XS{k}, 'LineWidth',1.5)
    title(compose('TH = %d', THv(k)));
    grid
end
function dXSdt=BSfunction(t,XS,TH)
X = XS(1); S = XS(2);
X0=2;
Sin=1000;
So=0;
S0=So+Sin;
V=20;
% TH=10;
ko=0.2;
kd=0.01;
Ks=150;
Y=0.5;
S=S0-(ko/Y)*X0*TH;
biomassgrowth=((ko*S*X)/(Ks+S))-(kd*X);
substratutilize = -((ko/Y)*X0);
dXSdt = [biomassgrowth; substratutilize];
end

The first figure plots them all on one axes, the second figure plots them each as separate subplots.

Best Answer

Related Solutions

MATLAB: System of arc-length defined ODEs with ode45

MATLAB: How to plot these 3 ode45 functions in 1 graph

Related Question