MATLAB Answers

euler method for solving system of ODE's 1st order

32 views (last 30 days)
Behzad Rahmani
Behzad Rahmani on 28 Feb 2021
Commented: Jan on 28 Feb 2021
Hello. I want to solve a system of first-order equations using the Euler method.
I have written this code but it does not give me the desired answer.Please check it and tell me the problem. thanks
a=0; %initial volume
b=10; %final volume
h = 0.1; % step
N = (b-a)./h;
nsteps = (b-a)/h + 1; % this is the number of elements in t(a:h:b) (It is = N+1)
f=zeros(1,nsteps);
g=zeros(1,nsteps);
v=a:h:b;
f = 10; % initial condition
g = 0;
p = 0;
q = 10;
F = @(v,f,g,p,q) -0.7*((200*(f/q))-((200^2)/50)*(g/q)*(p/q));
G = @(v,f,g,p,q) 0.7*((200*(f/q))-((200^2)/50)*(g/q)*(p/q));
P = @(v,f,g,p,q) -0.2*200*(p/q)+0.7*((200*(f/q))-((200^2)/50)*(g/q)*(p/q));
Q = @(v,f,g,p,q) f+g+p;
for i=1:N
f(i+1) = f(i) + h*F(v(i), f(i), g(i), p(i), q(i));
g(i+1) = g(i) + h*G(v(i), f(i), g(i), p(i), q(i));
p(i+1) = p(i) + h*P(v(i), f(i), p(i), p(i), q(i));
q(i+1) = q(i) + h*Q(v(i), f(i), q(i), p(i), q(i));
v(i+1)=a+i*h;
end
plot(v,f,'-',v,g,'--',v,p,'-o')
%plot(v,f,v,g);
  4 Comments
Jan
Jan on 28 Feb 2021
You still did not mention, why you assume that your code has a problem.

Sign in to comment.

Answers (1)

Alan Stevens
Alan Stevens on 28 Feb 2021
Simple Euler inaccurate for large step size. Reduce step size as in following and see if you get the output you expect:
a=0; %initial volume
b=10; %final volume
h = 0.005; % step Need a small step-size for simple Euler!!
N = (b-a)/h;
f = zeros(1,N+1); f(1) = 10;
g = zeros(1,N+1);
p = zeros(1,N+1);
q = zeros(1,N+1); q(1) = 10;
v=a:h:b;
F = @(f,g,p,q) -0.7*((200*(f/q))-((200^2)/50)*(g/q)*(p/q));
G = @(f,g,p,q) 0.7*((200*(f/q))-((200^2)/50)*(g/q)*(p/q));
P = @(f,g,p,q) -0.2*200*(p/q)+0.7*((200*(f/q))-((200^2)/50)*(g/q)*(p/q));
Q = @(f,g,p) f+g+p;
for i=1:N
f(i+1) = f(i) + h*F(f(i), g(i), p(i), q(i));
g(i+1) = g(i) + h*G(f(i), g(i), p(i), q(i));
p(i+1) = p(i) + h*P(f(i), p(i), p(i), q(i));
q(i+1) = q(i) + h*Q(f(i), q(i), p(i));
end
plot(v,f,'-',v,g,'--',v,p,'-o')
xlabel('v'), ylabel('f,g,p')
legend('f','g','p')
This results in:

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!