Fixing a plot to show a horizontal line instead of a point

Hello. First of all I don't know if my code is correct for the specific task, especially for the error calculations. My problem is that i get a point on the last graph on the errors figure but i need a horizontal line. Why does this happen ?

%--- Exercise 1 - A
%--- Solving the diff. eqn. N'(t) = N(t) - c*N^2(t) with c = 0.5 using
%--- Euler's method, Runge - Kutta method, and predictor-corrector method, and comparing in graph.
%--- I have four initial values for N0.

clc, clear all, close all

tmax = 100;
t = [0:tmax];
c = 0.5;
h = 0.1;
f = @(t,N) N - c*N.^2;
N0 = [0.1 0.5 1.0 2.0];
L = length(t)-1;
N = zeros(1,L);

for i = 1:length(N0)

%--- Exact solution ---

[T,n] = ode45(f,t,N0(i));

figure(1)
subplot(2,2,i), plot(T,n), hold on
title({'N'' = N - cN^2',sprintf('c = %.1f, N_0 = %.1f',c, N0(i))}),xlabel('t'), ylabel('N(t)'), hold on

%--- Euler's method ---

for j = 1:L
        
        NE(1) = N0(i);
        NE(1,j+1) = NE(j) + h*f(':',NE(j));
        ErrorEuler(j) = abs(n(j) - NE(j))/n(j)*100;
end

        plot(t,NE,'r-.'), hold on
        
%--- Runge - Kutta method ---
        
for j = 1:tmax
    
        NRK(1) = N0(i);
        
        K1 = f(t(j),NRK(j));
        K2 = f(t(j) + h*0.5, NRK(j) + h*K1*0.5);
        K3 = f(t(j) + h*0.5, NRK(j) + h*K2*0.5);
        K4 = f(t(j) + h, NRK(j) + h*K3);
    
        NRK(j+1) = NRK(j) + h*((K1 + K4)/6 + (K2 + K3)/3);
        
        ErrorRK(j) = abs(n(j) - NRK(j))/n(j)*100;
        
end

        plot(t,NRK,'b.'), legend({'Exact solution','Euler''s method','Runge - Kutta'},'Location','Southeast')
        
        figure(2)
        subplot(2,2,i)
        plot(NE(1:end-1),ErrorEuler,'.')
        hold on
        plot(NRK(1:end-1),ErrorRK,'-.')
        title({'N(t) vs. Error',sprintf('N_0 = %.1f',N0(i))}), xlabel('N(t)'), ylabel('Error (%)')
        legend({'Euler','Runge - Kutta'})
end

Answer 1

[Torsten]

N(t) = 2 for all t, and the error is 0 %. So plotting a line (vertical or horizontal) would be wrong in the case N_0 = 2 - you only get a single point.

[Left Terry]

Oh you are right, I though i was having time on the horizontal axis.Thanks!

[Torsten]

Your equation has an analytical solution. So instead of "ode45", I'd use "dsolve" to get a reference solution.

syms N0 c N(t)
eqn = diff(N,t) ==  N - c*N^2;
cond = N(0) == N0;
sol = dsolve(eqn,cond)

sol = 

fplot(subs(sol,[N0,c],[1,0.5]),[0 100])

●

[Left Terry]

@Torsten I know, I used dsolve, but when i used fplot to plot it, matlab couldn't plot it within the limits i gave. The graph was from 0 to 1 in both axis even when i set the correct limits. It's a bug or something, I couldn't fix it.

EDIT: It seems to have this problem when I put the value 0.5 in every form (0.5, 1/2, 5/10 etc) for c. I tried some other values and it works fine.

[Torsten]

In R2024b, it seems to work (see above).

Note that the value of h you set is wrong: h must be set to 1, not 0.1, for that the "exact" solution and your solutions (Euler/RK) can be compared.

[Left Terry]

@Torsten The problem states tthat we have to use h = 0.1.

[Walter Roberson]

The fplot() problem is specific to your release, R2016a. You can get around it by using xlim() and ylim()

[Left Terry]

@Torsten I tried it and it didn't work.

[Torsten]

@Left Terry The problem states tthat we have to use h = 0.1.

But you use h = 1 in the ode45 calculation:

tmax = 100;
t = [0:tmax]

t = 1×101

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

So either you use h = 0.1 or h = 1 in both calculations.