eulers method and relative error

hi im new to matlab and have this question ,

Use Euler’s method to evaluate the solution to the system of differential equations over the first 20 hours. Use a step size of h = 15 minutes and examine the relative error between the numerical and analytical calculations.

these are the equations

1) R + k*x2(t) - k*x1(t)

2) k*x1(t) - (alpha + k) * x2(t)

values of aplha = 10

k = 50

and R = 100

can anyone help ?

[Jim Riggs]

What you need to do is write a program that implements Euler's method to solve both of these equations, and compare the result to the analytical (i.e. exact) answer. This will give you an measure of the amount of numerical error in the Euler's method.

[courteney fishwick]

hi thanks ,this is the code ive wrote but getting errors .

Array indices must be positive integers or logical values.

Error in sym/subsref (line 902)

R_tilde = builtin('subsref',L_tilde,Idx);

Error in ppoorrtteeuull (line 26)

ics = [X1_0 == X1(0); X2_0 == X2(0)];

% Numerical solutions to ODEs: Eulers method [part1]

clear all 

syms f1(t,X1,X2) f2t(t,X1, X2) x1(t) x2(t) tpt x1pt x2pt Dx1pt Dx2pt rel_error k  R
R = 100
k = 50
aplha = 10 
% RHS of general ODE: Dx(t) = f1(t,X,Y) , DY(t) = f2(t,X,Y) and initial condition [part2] 
f1(t,X1,X2) = R + k*x2(t) - k*x1(t)
f2(t,X1,X2) = k*x1(t) - (aplha + k) * x2(t)
X1_0 = 150
X2_0 = 150

% Eulers method : h=step size and N=number of iterations
h = 0.25
N = 20

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Differential equation 
deq1 = diff( x1(t) , t ) == f1(t, X1, X2)
deq2 = diff( x2(t) , t ) == f2( t, X1, X2)

deq = [deq1; deq2];
ics = [X1_0 == X1(0); X2_0 == X2(0)];

SOL = dsolve(deq, ics)
X1sol = @(T) subs(SOL.x1(t), t, T);
X2sol = @(T) subs(SOL.x2(t), t, T);

for k=1:N+1
    tpt(k) = h*(k-1);
    if k==10
        x1pt(k) = x150; % Initial condition
        x2pt(k) = y150; % initial condition 
    else 
        x1pt(k) = x1pt(k-1) + h*Dx1pt(k-1); % Eulers formula
        x2pt(k) = x2pt(k-1) + h*Dx2pt(k-1); % Eulers formula
    end
    Dxpt(k) = f1( tpt(k), x1pt(k), x2pt(k) ); % f1(t,X,Y)
    Dypt(k) = f2( tpt(k) ,  x1pt(k), x2pt(k) ); % f2(t,X,Y)
    X1rel_error(k) = 100 * abs( ( X1sol(tpt(k)) - x1pt(k) ) /X1sol(tpt(k)) ); % Relative error for xpt
    X2rel_error(k) = 100 * abs( ( X1sol(tpt(k)) - x2pt(k) ) /X2sol(tpt(k)) ); % Relative error for ypt
end

X1REL_ERROR = vpa(xrel_error, 10);
x2REL_ERROR = vpa(yrel_error, 10);

TMAX = h*N;

subplot(2,2,1)
fplot( @(t) X1sol(t), [0 TMAX],...
    'color' , 'r',...
    'linewidth' , 1.3)
grid on
title('TBC')
xlabel('time, t ')
ylabel('TBC')

hold on 
plot( tpt, xpt,...
    'color', 'b',...
    'linestyle', '-',...
    'marker', 'o',...
    'markersize', 6,...
    'markerfacecolor', 'k',...
    'markeredgecolor', 'k')
grid on 
legend('TBC')
hold off

[Jim Riggs]

Well, I don't use the symbolic processing toolbox, so I don't have any advice regarding that.

But the reference to X1(0) and X2(0) may be the problem. If the numbers in the parenteses are subscripts, they must be positive, that would be what the error is indicating. A value of zero is not positive. Needs to be 1 or greater (integer or logical type)