script Runge Kutta4

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

//Programa para resolver pelo método de Euler o Sitema: x1'(t)+x2(t)=t^2, x2'(t)-x1(t)=1 com 
//as condições iniciais x1(0)=1, x2(0)=0
t0=0 //tempo inicial
tf=2*%pi //tempo final
h=sqrt(sqrt(150*10^(-5)/(exp(tf-t0)-1)))//passo necessário para que o erro global seja inferior a 5e-5
//vou definir a matriz X=(x(i,j)) onde i=1,2 e j=1,...,n (n= numero de iterações)
//note que a primeira iteração é indexada com 1 não com zero
//h=5*h //só pra diminuir o tamanhao da matriz X
X(1,1)=1; //valor exato de x1 na 1º aproximação
X(2,1)=0; //valor exato de x2 na 1º aproximação
i=1
t=t0
Errox1=0
Errox2=0
Auxiliar=0
//Y(1,1)=1;
//Y(2,1)=0;
while t<tf
//vetor m1
m11=(-X(2,i)+t^2);
m12=(1+X(1,i));
//vetor m2
m21=(-X(2,i)-h*m12/2+(t+h/2)^2);
m22=(1+X(1,i)+h*m11/2);
//vetor m3
m31=-(X(2,i)+h*m22/2)+(t+h/2)^2;
m32=X(1,i)+h*m21/2+1;
//vetor m4
m41=-(X(2,i)+h*m32)+(t+h)^2;
m42=X(1,i)+h*m31+1;
X(1,i+1)=X(1,i)+h*(m11+2*m21+2*m31+m41)/6;
X(2,i+1)=X(2,i)+h*(m12+2*m22+2*m32+m42)/6;
Errox1=(-2*sin(t+h)+2*cos(t+h)+2*(t+h)-1)-X(1,i+1);
Errox2=(2*sin(t+h)+2*cos(t+h)+(t+h)^2-2)-X(2,i+1);
 if abs(Errox1) > abs(Errox2) then
 Emax=Errox1;
 elseif abs(Errox2) > abs(Errox1)
 Emax=Errox2;
 else
 Emax=Errox2;
 end
 if abs(Emax)<abs(Auxiliar) then
 Emax=Auxiliar;
 else
 Emax=Emax;
 end
Auxiliar=Emax;
t=t+h;
i=i+1;
end
T=[0:h:2*%pi+h];
A=X';
X1=A(:,1); X2=A(:,2);
ExataX1=[-2*sin(T)+2*cos(T)+2*T-1]';
ExataX2=[2*sin(T)+2*cos(T)+T^2-2]';
X1eExata=[X1, ExataX1];
X2eExata=[X2, ExataX2];
xsetech([0,0,0.5,0.5]); xtitle("X1 aproximado e exato"); plot2d(T,X1eExata, style=[-1,1],...
axesflag=5,leg="x1 aproximado@x1 exato")
xsetech([0.6,0,0.5,0.5]); xtitle("X2 aproximado e exato"); plot2d(T,X2eExata, style=[-1,1], ...
axesflag=5,leg="x2 aproximado@x2 exato")