euler method for solving system of ODE's 1st order
35 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (1)
Alan Stevens
am 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 Kommentare
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!