Skip to main content

SIMPLE AND MODIFIED EULER METHOD


##Usage:Call Octave from terminal
##and then call EulerMethodUmitAlkus.m
##from octave and finally
##press enter. That's all.
##Simple Euler Method
##Constants and initializations
x=[]; ## initial empty vector for x
y=[]; ## initial empty vector for y
x(1)=1; ## initial value of x
y(1)=1; ## initial value of y
h=1E-3; ## increment in x
dery=[]; ## 1st derivative of y wrt x
dery(1)=0;## 1st entry of dery
n=1; ## inital loop index for while
## enter the while loop for the interval x=[1,2]
while (x(n)<=2)
x(n+1)=x(n)+h;
dery(n+1)=x(n)*x(n)-2*y(n)/x(n); ##given
y(n+1)=y(n)+h*dery(n+1); ##Euler method
n++;
endwhile
##exit from the 1st while loop
##Modified Euler Method
##Constant and initializations
x(1)=1; ## beginnig of the interval [1,2]
ymod(1)=1; ## inital value for modified y.
ymid=[]; ## empty vector function evaluated at x midpoint
xmid=[]; ## empty vector func. of midpoints of the interval h in x-axis.
derymod=[]; ## modified derivatives of ymod wrt x.
derymid=[]; ## derivative of ymid wrt xmid.
derymid(1)=0; ##1st entry of der. of y wrt x at midpoints.
##Enter the 2nd while loop
n=1;
while (x(n)<=2) ## x(n)'s are the beginning values of the interval h.
xmid(n)=x(n)+h/2;
derymod(n)=x(n)*x(n)-2*ymod(n)/x(n); ##given equation
ymid(n)=ymod(n)+derymod(n)*h/2;
derymid(n+1)=xmid(n)*xmid(n)-2*ymid(n)/xmid(n);
ymod(n+1)=ymod(n)+derymid(n)*h; ##modified Euler Method
x(n+1)=x(n)+h;
n++;
endwhile
## Plot for Simple Euler Method
subplot(2,1,1)
hold on
plot(x,dery,'c-');
plot(x,ymod,'r-');
legend('dy/dx','y')
title('Simple Euler Method');
xlabel('x');
hold off
##Plot for Modified Euler Method
subplot(2,1,2)
hold on
plot(x,derymid,'m-');
plot(x,ymod,'g-');
title('Modified Euler Method');
xlabel('x');
legend('dy/dx(mid)','y-mod');
hold off
print('-dpsc','SIMMODEULER.ps')
save -text SIMMODEULER.dat

Comments

Popular posts from this blog

NEWTON-RAPSON METHOD FOR HEAT FLOW

##Constants and initializations
a=5.67E-8; ## Stefan-Boltzman constant[Watt/meter^2Kelvin^4]
e=0.8; ## Rod surface emissivity [Dimensionless]
h=20; ## Heat transfer coefficient of air flow [W/m^2-K]
Tinf=Ts=25; ## Temperature of air and the walls of the close[Celcius]
D=0.1; ## Diameter of the rod[meter]
I2R=100; ## Electric power dissipated in rod (Ohmic Heat)[W]
T=[]; ## Temperature of the rod[*C]
T(1)=25; ## Initial guess of the temperature of the rod[*C]
Q=[]; ## Heat function [W]
Qp=[]; ## First derivative of Q wrt T [W/C*].
for i=1:100
Q(i)=pi*D*(h*(T(i)-Tinf)+e*a*(T(i)^4-Ts^4))-I2R;
Qp(i)=pi*D*(h+4*e*a*T(i)^3);
T(i+1)=T(i)-Q(i)/Qp(i); ## Newton-Rapson Method
endfor
printf('The steady state temperature is %f\n',T(i+1))
save -text HeatFlowTemp.dat
## The plot
t=1:100; ##temperature
for n=1:100
H(n)=pi*D*(h*(t(n)-Tinf)+e*a*(t(n)^4-Ts^4))-I2R;
endfor
plot(t,H)
xlabel('T(Celcius)');
ylabel('Q(Watt)');
legend('Q(T)');
title('Heat flow vs Temperatu…

INITIAL POSITIONS OF DIFFUSIVE PARTICLES

## A function takes in the number of diffusive particle
## as an argument and gets the initial positions in 2D
## as an outcome by inserting them in a square.
function p=initial_positions7(Nwalkers)
dt=3; ## to extend the interval of axes.
a=sqrt(Nwalkers); ## step size.
n=(a-1)/2; ## max range in the x and y axes.
if(rem(a,2)~=1) ## Nwalkers must be odd squared.
printf("The number you have entered is not an odd squared. Exiting...\n");
return
endif
n=(a-1)/2; ## max range in the x and y axes.
x=[-n:n]'; ## integer interval in the x axis.
p=zeros(a,2); ## initialize the positions of 'a' walkers.
p(:,1)=x; ## equate the column of p to x vector.
p=repmat(p,a,1); ## replicate positions of 'a' walkers to get a*a walkers.
for i=1:a:Nwalkers ## the loop through the number walkers.
m=(i+a-1)/a; ## index get the entries of x resp.
p(i:i+a-1,2)=x(m) ; ## y components are a-fold degenerate.
endfor
plot(p(:,1),p(:,2),'b*;;')
axis(dt*[-n, n,-n,n])
endfunction