Profile icon
a5rocks

Subject yourself to a pretty terrible great programming language. Uncomment "#example.o", and edit the file "example.o" to make your own scripts.

You are viewing a single comment. View All
Profile icon
GabriellaNeves

function principal
t0 = 0;
tf = 50;
x0 = 99000;
k1 = 2*(10^(-6));
k2 = 10^(-4);
m = 100000;
options = odeset('Abstol', k1,'Reltol', k2);
[t y] = ode45(@fdo, [t0 tf], x0, options);
plot(t,y)
function dz = fdo(t,z)
k1 = 2*(10^(-6));
k2 = 10^(-4);
x0 = 99000;
m = 100000;
dz = k2*(m - z - x0 * exp(-(k1/k2)z));
% function y = euler(fun,alpha,a,b,N,grau)
%
% y = euler(fun,alpha,a,b,N,grau)
%
% METODO DE EULER: resolve um EDO de primeira ordem
% com uma condicao inicial y(a)=alpha.
%
% fun - funcao continuamente diferenciavel
% alpha - condicao inicial
% [a,b] - intervalo
% N - numero de passos
% Default h=0.01
% grau - grau do metodo Euler, grau pode ser 1, 2 ou 3.
% Default grau=1
%
% quero x(t), y(t) e z(t)
syms z
m = 100000;
x = 99000;
k1 = 2 * 10^(-6);
k2 = 10^(-4);
t = [0:50];
h = 0.01;
%dz = k2
(m - z - x0 * exp(-(k1/k2)z));
%x = x0
exp^(-(k1/k2)z);
%y = m - x - z;
N = 12.5;
alpha = 99000;
grau = 4;
y = alpha;
fun = k2
(m - z(t) - x*(exp^(-(k1/k2)*z(t))));
if grau==1

for k = 2:(N+1)
phi = feval(fun,t,x,0);
t = t + h;
y(t) = t;
x = x + h*phi;
x(t) = x;
end
end
if grau==2

for k = 2:(N+1)
phi = feval(fun,t,x,0) + (h/2)[feval(fun,t,x,1)];
t = t + h;
y(t) = t;
x = x + h
phi;
x(t) = x;
end
end
if grau==3

for k = 2:(N+1)
phi = feval(fun,t,x,0) + h/2*[feval(fun,t,x,1)] +
(h^2)/6*[feval(fun,t,x,2)];
t = t + h;
y(t) = t;
x = x + h*phi;
x(t) = x;
end
end
if grau>=5