fork download 
  1. % Define the differential equation dy/dx = f(x, y)
  2. %euler loop
  3. f = @(x, y) 2*x;   %  f(x, y) = 2*x
  4.  
  5. % Step size and number of steps
  6. h = 0.4;                   % step size
  7. n = (2 - 0)/h;             % total number of steps to reach x = 2 from x = 0
  8.  
  9. % Set initial values
  10. x(1) = 0;                  % initial x0
  11. y(1) = 1;                  % initial y0
  12.  
  13. % Euler method loop
  14. for i = 1:n
  15.     x(i+1) = x(i) + h;             % update x to next step
  16.     y(i+1) = y(i) + h * f(x(i), y(i));  % Euler formula: y_{i+1} = y_i + h*f(x_i, y_i)
  17. end
  18.  
  19. %modified euler method
  20.  
  21. %write the function
  22.  
  23. f = @(x, y) exp(x^2);
  24.  
  25. %x,y which makes two inputs, here y is not use 
  26.  
  27. %step size
  28.  
  29. h = 0.1;
  30. n = (2-0)/h;
  31. %initial value
  32. x(1) = 0;
  33. y(1) = 1;
  34.  
  35. for i = 1:n
  36. x(i+1) = x(i) + h; %step size
  37. %euler method
  38. y(i+1) = y(i) + h*f(x(i), y(i));
  39. %mod euler
  40. y(i+1) = y(i) + h * (f(x(i), y(i)) + f(x(i+1), y(i+1))) / 2;
  41.  
  42. end 
  43.  disp([x' y'])
  44.  
  45.  
  46. %  e^x^2 from 0 to 1 simpson 1/3 n=12(n should be even)
  47.  
  48. f=@(x) exp(x.^2);
  49.  
  50. a=0;
  51. b=1;
  52.  
  53. n=12;
  54.  
  55. h=(b-a)/n;
  56.  
  57. sum=f(a)+f(b);
  58.  
  59. %loop niboh
  60. for i=1:n-1
  61.     a=a+h;
  62.     if mod(i,2)==1
  63. %odd hobe 
  64.         sum=sum+4*f(a); 
  65.     else
  66. %even hobe bhag kore jdi dkhi
  67.         sum=sum+2*f(a);
  68.     end
  69. end
  70. Sum=h/3*sum
  71.  
  72. fprintf('Approximate integral = %.10f\n',Sum);
  73.  
  74.  
  75. %simpson 3/8 
  76. clear all;
  77. clc;
  78. format long;
  79. %function boshate hobe
  80.  
  81. f=@(x) exp(x.^2);
  82.  
  83. %initial value boshate hobe
  84. a=0;
  85. b=1;
  86.  
  87. n=12;
  88.  
  89. %step size ber korbo
  90. h=(a-b)/n;
  91.  
  92. %formula onujai x0 are xn jug hobe total with last and first value
  93.  
  94. sum= f(a)+f(b);
  95.  
  96. %loop dibo
  97. for i=1:n-1
  98.   a=a+h; %x er next value ber korbo
  99.  
  100.   if mod(i,3)==0 %i multiple 3 er 
  101.      sum=sum +2*f(a); %3 er gunitok gula ashbe
  102.  
  103.      else
  104.      sum=sum+3*f(a);
  105.      end 
  106.  
  107.     end  
  108.  
  109.     %ekhn just last er 3h/8 diye kaj shesh krbe
  110.     sum=3*h*sum/8;
  111.     fprintf('Approximate integral = %.10f\n', sum);
  112.  
  113.  
  114. %TROPZOIDAL e^x^2 from tropzoidal n=12
  115.  
  116. %function likhbo
  117. f=@(x) exp(x.^2);
  118.  
  119. %upper lowe limit
  120. a=0;
  121. b=1;
  122.  
  123. %n er value likhbo
  124. n=12
  125.  
  126. %step size ber krbo
  127.  
  128. h=(b-a)/n;
  129.  
  130. %sum krbo first and last
  131. sum=f(a)+f(b);
  132.  
  133. %loop dibo
  134. for i=1:n-1
  135.   %x er next value ber krbo
  136.   a=a+h;
  137.   sum=sum+2*f(a);
  138.   end 
  139.   sum=h/2*sum
  140.   fprintf('Approximate intregal: %.10f\n', sum')
  141.  
  142.   %for exact value ber korbo abar function likhbo
  143.   f= @(x) exp(x.^2);
  144.   result=quad(f,0,1)
  145.  
  146.  
  147. %midpoint
  148.  
  149. %dy/dx=x+y y(0)=1 where step size h=0.1 
  150. %function likhbo
  151. f=@(x,y) x+y;
  152. %from question x=0 y=1 
  153. x(1)=0;
  154. y(1)=1;
  155. % x er man dewa 0.2
  156. xn=0.2;
  157. h=0.1;
  158.  
  159. n=(xn-x(1))/h;
  160.  
  161. %loop dibo
  162. for i=1:n 
  163.  
  164. %midpoint formula dibo
  165. k1=f(x(i),y(i));
  166. k2=f(x(i)+h/2,y(i)+h/2*k1);
  167. y(i+1)=y(i)+h*k2;
  168. x(i+1)=x(i)+h;
  169.  
  170. end 
  171. [x' y']
  172.  
  173.  
  174. %ratson method
  175.  
  176. %dy/dx=x+y y(0)=1 where step size h=0.1 
  177. %function likhbo
  178. f=@(x,y) x+y;
  179. %from question x=0 y=1 
  180. x(1)=0;
  181. y(1)=1;
  182. % x er man dewa 0.2
  183. xn=0.2;
  184. h=0.1;
  185.  
  186. n=(xn-x(1))/h;
  187.  
  188. %loop dibo
  189. for i=1:n 
  190.  
  191. %ralston formula dibo
  192. k1=f(x(i),y(i));
  193. k2=f(x(i)+3*h/4,y(i)+3*h*k1/4);
  194.     y(i+1)=y(i)+h*(k1+2*k2)/3;
  195.     x(i+1)=x(i)+h;
  196.  
  197. end 
  198. [x' y']
  199.  
  200.  
  201. %rk 4th method
  202.  
  203. %dy/dx=x+y y(0)=1 where step size h=0.1 
  204. %function likhbo
  205. f=@(x,y) x+y;
  206. %from question x=0 y=1 
  207. x(1)=0;
  208. y(1)=1;
  209. % x er man dewa 0.2
  210. xn=0.2;
  211. h=0.1;
  212.  
  213. n=(xn-x(1))/h;
  214.  
  215. %loop dibo
  216. for i=1:n 
  217.  
  218.    % R-K 4th order
  219.      k1=f(x(i),y(i));
  220.      k2=f(x(i)+h/2,y(i)+h*k1/2);
  221.     k3=f(x(i)+h/2,y(i)+h*k2/2);
  222.     k4=f(x(i)+h,y(i)+h*k3);
  223.     y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6;
  224.     x(i+1)=x(i)+h;
  225.  
  226. end 
  227. [x' y']
  228.  
  229.  
  230. %heuns
  231.  
  232. %dy/dx=x+y y(0)=1 where step size h=0.1 
  233. %function likhbo
  234. f=@(x,y) x+y;
  235. %from question x=0 y=1 
  236. x(1)=0;
  237. y(1)=1;
  238. % x er man dewa 0.2
  239. xn=0.2;
  240. h=0.1;
  241.  
  242. n=(xn-x(1))/h;
  243.  
  244. %loop dibo
  245. for i=1:n 
  246.  
  247.    %  Heun's
  248.     %k1=f(x(i),y(i));
  249.     %k2=f(x(i)+h,y(i)+h*k1);
  250.     %y(i+1)=y(i)+h*(k1+k2)/2;
  251.     %x(i+1)=x(i)+h;
  252.  
  253. end 
  254. [x' y']
  255.  
  256.  
  257. %graph of euler loop
  258. % Define the function dy/dx =2x= f(x, y)
  259. f = @(x, y) exp(x^2);
  260.  
  261. % Step size and number of steps
  262. h = 0.001;
  263. xn=2;
  264. n = (2 - 0) / h;
  265.  
  266. % Set initial values
  267. x(1) = 0;
  268. y(1) = 1;
  269.  
  270. % Euler loop
  271. for i = 1:n
  272.     x(i+1) = x(i) + h;
  273.     y(i+1) = y(i) + h * f(x(i), y(i));
  274.     y(i+1)=y(i)+h*(f(x(i),y(i))+f(x(i+1),y(i+1)))/2;  %modi-Euler
  275. end
  276. n
  277. % Exact solution
  278. y_exact = x.^2 + 1;
  279. % Plotting both
  280. plot(x, y, 'ro--', x, y_exact, 'b.--');
  281. legend('Euler Approx', 'Exact Solution');
  282.  
  283.  
  284. % weddles n=12 (n shold be multiple of 6)
  285.  
  286. f = @(x) exp(x.^2);
  287.  
  288. a = 0;
  289. b = 1;
  290. n = 12; % subintervals, must be multiple of 6 for Weddle's rule
  291.  
  292. h = (b - a) / n;
  293.  
  294. % Generate points
  295. x = a:h:b;
  296. y = f(x);
  297.  
  298. sum = 0;
  299.  
  300. % Apply Weddle's rule in blocks of 6 subintervals (7 points)
  301. for k = 1:6:n
  302.     sum = sum + (3*h/10) * ( ...
  303.         y(k)   + 5*y(k+1) + y(k+2) + ...
  304.         6*y(k+3) + y(k+4) + 5*y(k+5) + y(k+6) );
  305. end
  306.  
  307. sum
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
Success #stdin #stdout 0.54s 57728KB
comments (?)
stdin 
Standard input is empty
stdout 
    0.00000    1.00000
    0.10000    1.10050
    0.20000    1.20305
    0.30000    1.30979
    0.40000    1.42318
    0.50000    1.54606
    0.60000    1.68192
    0.70000    1.83521
    0.80000    2.01165
    0.90000    2.21887
    1.00000    2.46717
    1.10000    2.77076
    1.20000    3.14947
    1.30000    3.63148
    1.40000    4.25742
    1.50000    5.08677
    1.60000    6.20795
    1.70000    7.75441
    1.80000    9.93076
    1.90000   13.05575
    2.00000   17.63396
Sum =  1.4627
Approximate integral = 1.4626661264
Approximate integral = -1.4626836999
n =  12
sum =  1.465794273401113
Approximate intregal: 1.4657942734
result =  1.462651745907181
ans =

   0.000000000000000e+00   1.000000000000000e+00
   1.000000000000000e-01   1.110000000000000e+00
   2.000000000000000e-01   1.242050000000000e+00

ans =

   0.000000000000000e+00   1.000000000000000e+00
   1.000000000000000e-01   1.110000000000000e+00
   2.000000000000000e-01   1.242050000000000e+00

ans =

   0.000000000000000e+00   1.000000000000000e+00
   1.000000000000000e-01   1.110341666666667e+00
   2.000000000000000e-01   1.242805141701389e+00

ans =

   0.000000000000000e+00   1.000000000000000e+00
   1.000000000000000e-01   1.110341666666667e+00
   2.000000000000000e-01   1.242805141701389e+00

n =  2000
sum =  1.462652067662139
https://ideone.com/7JTyvb
language:
Octave (octave 4.4.1)
created:
2 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!