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

FACTORIAL

## Function that calculates the factorial of a number ## Usage : f=factorial(n) function f=factorial(n) ## Initialize the output f=1; ## Check whether the input is correct if ( (n<0) || (rem(n,1)~=0) ) printf("n cannot be a negative number. Exiting...\n"); return endif for num=1:n f*=num; endfor endfunction

One Dimensional Harmonic Oscillator-Numerov Method

x=[]; h0=1; M=4; N=M+1; x(1)=0; x(N)=x(1)+h0*M; x=x(1):h0:x(N) A=zeros(N); A(1,1)=-2*(5*(x(N)*h0)^2/12+1); A(N,N)=-2*(5*(x(1)*h0)^2/12+1); A(1,2)=1-(x(M)*h0)^2)/12; A(N,M)=1-(x(2)*h0)^2)/12; B=zeros(N); B(1,1)=B(N,N)=-10*(h0^2)/6; B(1,2)=B(N,M)=-(h0^2)/6; for i=2:M B(i,i)=-10*(h0^2)/6; B(i,i-1)=B(i,i+1)=-(h0^2)/6; A(i,i)=-2*(5*(x(N+1-i)*h0)^2/12+1); A(i,i+1)=1-(x(N-i)*h0)^2)/12; A(i,i-1)=1-(x(N+2-i)*h0)^2)/12; end A B

REALISTIC AND NON-REALISTIC STRINGS

## A 'realistic', 'non-elastic' string, which responses to any ## bending and has stifness. This script takes in the previous ## and the present profiles and iterates to find ## the profile in the next time step. The ratio 'r' is not 1 ## like in the 'non-realistic' string since the speed of the wave ## always less than the speed of the string it should be less than 1 ## for best and most stable solution ##constants dx=1e-2 ## Spatial increment (m) L=2 ## Length of the string (m) M=L/dx ## Dimensionless partition E=1e-4 ## Dimensionless stiffnes function ynext=propagate_stiff(ynow,yprev,r) ## Quick and dirty way to fix boundary conditions -- for each step ## they are the same as the previous step. ynext=ynow; ynow(1)=ynow(2)=0; ##Entering the loop for i=3:length(ynow)-1 ## boundary condition ynow(length(ynow)-1)=ynow(length(ynow)-2)=0; ## Divide the ynext with many terms into three parts for easiness ynext(i)=(2−(2*r^2)−(6*E*(r^2)*...