Runge-kutta third order method:
%rk3:runge kutta of thirdorder
clc;
clear all;
close all;
% y' = y-x ode condition
f = @(x,y) y-x;
fex = @(x) exp(x)+x+1; % exact solution
a=0;
b= 3.2;
n =16;
h=(b-a)/n;
y(1) =2; %initial value
i = 0;
for x= a:h:b
i = i+1;
K1 = f(x,y(i)); %initializing solution
K2 = f(x+h*0.5,y(i)+h*K1*0.5);
K3 = f(x+h, y(i)-h*K1 +2*K2*h);
y(i+1) =y(i)+h*(1/6)*(K1 +4*K2+K3);
g(i) = fex(x);
xx(i) = x;
Error(i) = abs(g(i) - y(i)); %error obtain
end
%plot result
plot(xx,y(1:n+1),'k',xx,g,'y')
legend('RK3','Exact solution')
xlabel('x')
ylabel('y')
title('RK3 vs exact solution')
I am getting wrong error value.. please check my code

4 commentaires

KSSV
KSSV le 15 Mai 2019
Wrong error value????
NOTE: initialze the variable inside the loop.
SHIVANI TIWARI
SHIVANI TIWARI le 15 Mai 2019
am not getting the right value of absolute error.
i have taken same problem 2. but am not getting the right value of absolute error.
and plz also my program for butcher 5 for solving 1st order ode
clc;
clear all;
close all;
% y' = y-x ode condition
f = @(x,y) y-x;
fex = @(x) exp(x)+x+1; % exact solution
a=0;
b= 3.2;
n =16;
h=(b-a)/n;
y(1) =2; %initial value
i = 0;
for x= a:h:b
i = i+1;
K1 = f(x,y(i));
K2 = f(x+h*(1/4),y(i)+h*(1/4)*K1);
K3 = f(x+h*0.5, y(i)+ h*K1*(1/8)+h*K2*(1/8));
K4 = f(x+h*0.5,y(i)- h*K2*0.5+h*K3);
K5 = f(x+3*h*(1/4),y(i)+ (3/16)*h*K1+(9/16)*h*K4);
K6 = f(x+h,y(i)- 3*h*K1*(1/7)+2*h*K2*(1/7) + (12/7)*h*K3 -(12/7)*h*K4 + (8/7)*h*K5);
y(i+1) =y(i)+h*(1/90)*(7*K1 +32*K3+12*K4+ 32*K5 +7*K6);
g(i) = fex(x);
xx(i) = x;
Error(i) = abs(g(i) - y(i)); %error obtain
%plot result
plot(xx,y(1:n+1),'r',xx,g,'g','LineWidth',2)
legend('butcher5','Exact solution')
xlabel('x')
ylabel('y')
title('butcher5 vs exact solution')
excuse me, how do you determine those formulas for finding the values for K's?
SHIVAM
SHIVAM le 6 Nov 2023
thank u so much

Connectez-vous pour commenter.

Réponses (0)

Catégories

En savoir plus sur Numerical Integration and Differential Equations dans Centre d'aide et File Exchange

Commenté :

le 6 Nov 2023

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by