euler method for solving system of ODE's 1st order
34 views (last 30 days)
Show older comments
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
Answers (1)
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:
0 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!