1. POLYFIT
x=[ 10 12 16 20 25 30];
y=[29 33 41 53 65 70];
p=polyfit(x,y,1);
x1=linspace(x(1),x(end));
y1=polyval(p,x1);
plot(x,y,'*',x1,y1)
legend('data points','linear curve')
2. Dsolve
y=dsolve('Dy=3*x+y/2','y(0)=1','x')
x=0:0.1:1;
z=eval(y);
plot(x,z,'*')
2.a ODE
y0=1;
xspan=[0,1];
f=@(x,y) 3*x+y/2;
[x,y]=ode45(f,xspan,y0);
plot(x,y,'r')
2: Euler’s Method
function [x,y] = EulerODE(x0, xn, y0, h, f)
x=x0:h:xn;
n=length(x)-1;
y=zeros(1,n+1);
y(1) = y0;
for j=1:n
y(j+1) = y(j) + h*f(x(j), y(j));
end
end
3. Euler’s Modified Method
function []= modEulerODE(x1, y1, h)
f=@(x,y) 3.*x+y/2; %put function
xn=0.2; %change value here
x = x1: h: xn;
n=length(x);
y(1)=y1;
maxit = 10;
for i=2:n
yp(i)=y(i-1)+h*f(x(i-1), y(i - 1));
s(i, 1)=yp(i);
for j=2:maxit
s(i,j)=y(i-1)+(h/2)*(f(x(i-1), y(i-1)) + f(x(i), s(i, j-1)));
if abs(s(i, j)-s(i,j-1)) <10^(-6)
y(i)=s(i,j);
break
end
end
end
fprintf('Approximate. sol. y(%f)=%f', x(end),y(end))
plot(x, yp, 'r--')
hold on
plot(x, y, 'b*')
hold on
plot(x,y,'r--')
hold on
u=dsolve('Dy=3*x+y/2', 'y(0)=1', 'x');
u=eval(u);
plot(x, u, 'g--')
legend('Approximate. sol. by Euler', 'Appr. sol. by Modified Euler', 'Exact sol.')
end
4. Runge-Kutta Fourth Order Method
function [x,y]=RK4ODE(x0,xn, y0,h,f )
x = x0:h:xn;
n = length(x)-1;
y = zeros(1,n+1);
y(1) = y0;
for i=1:n
K1 = h*f(x(i),y(i));
K2 = h*f(x(i)+0.5*h,y(i)+0.5*K1);
K3 = h*f(x(i)+0.5*h,y(i)+0.5*K2);
K4 = h*f(x(i)+h,y(i)+K3);
y(i+1) = y(i)+(1/6)*(K1+2*K2+2*K3+K4);
end
end
5. Milne’s Method
function [x,y]=MilneODE(x0,xn, y0,h,f )
x=x0:h:xn;
n=length(x)-1;
y = zeros(1,n+1);
y(1)=y0;
for i=1:3
K1=h*f(x(i),y(i));
K2=h*f(x(i)+0.5*h,y(i)+0.5*K1);
K3=h*f(x(i)+0.5*h,y(i)+0.5*K2);
K4=h*f(x(i)+h,y(i)+K3);
y(i+1)=y(i)+(1/6)*(K1+2*K2+2*K3+K4);
end
for i=4:n
y(i+1)=y(i-3)+(4*h/3)*(2*f(x(i),y(i)) - f(x(i-1),y(i-1))+2*f(x(i-2),y(i-2))); %predictor formula
y(i+1)=y(i-1)+(h/3)*(f(x(i+1),y(i+1)) + 4*f(x(i),y(i)) + f(x(i-1),y(i-1))); % Corrector formula
end
6. Script files to run any of the above functions
Note: Keep the script file and function file in the same folder.
x0=0; xn = 1; y0=0.5; h=0.1; f = @(x,y) y-x^2+1; % Change according to the problem.
% Uncomment any one of the below line to the output for the given problem.
%[x,y] = EulerODE(x0, xn, y0, h, f);
[x,y] = modEulerODE(x0, xn, y0, h, f);
%[x,y] = RK4ODE(x0, xn, y0, h, f);
%[x,y] = MilneODE(x0, xn, y0, h, f);
% Display the values of x and y using 'disp' or 'fprintf'
fprintf('Approximate solution of the ODE is\n');
disp([x;y]);
ye=eval(dsolve('Dy=y-x^2+1','y(0)=0.5','x'));% Change according to the problem
plot(x,y,'*',x,ye,'r');
Run the script file ( Saved as ODEScript.)
>> ODEScript
Approximate solution of the ODE is
0 0.2500 0.5000 0.7500 1.0000
0.5000 0.9205 1.4256 2.0039 2.6408
7. One Dimensional Heat Equation (Explicit method)
function oneD_heat(t0, tn,x0,xn,h,k,c )
t=t0:k:tn;
x=x0:h:xn;
m=length(x);
n=length(t);
a=c*k/h^2
f=@(x) sin(pi*x); % Change f according with the problem
u=zeros(m,n);
u(:,1)=f(x);
if a> 0.5
fprintf('The method fails')
return
end
for j=1:n-1
for i=2:m-1
u(i, j+1)=a*u(i-1,j)+(1-2*a)*u(i,j)+a*u(i+1,j);
end
end
for j=1:n
plot(x,u(:,j))
hold on
end
figure
surf(t,x,u)
xlabel('t')
ylabel('x')
zlabel('u')
end
8.Adams Bash forth Method
function [x,y]=Adams_method(x0,y0,xn,n)
h=(xn-x0)/n;
f=@(x,y) 0.09*y-20000;
x=x0:h:xn;
y=zeros(1,n);
maxit=10;
s(:,:)=zeros(n,maxit);
y(1)=y0;
for i=1:3
K1=h*f(x(i),y(i));
K2=h*f(x(i)+0.5*h,y(i)+0.5*K1);
K3=h*f(x(i)+0.5*h,y(i)+0.5*K2);
K4=h*f(x(i)+h,y(i)+K3);
y(i+1)=y(i)+(1/6)*(K1+2*K2+2*K3+K4);
end
%predictor formula
for i=4:n
y(i+1)=y(i)+(h/24)*(55*f(x(i),y(i))-...
59*f(x(i-1),y(i-1))+37*f(x(i-2),y(i-2))-...
9*f(x(i-3),y(i-3)));
s(i+1,1)=y(i+1);
%corrector formula
for j=1:maxit
s(i+1,j+1)=y(i)+(h/24)*(f(x(i-2),y(i-2))-...
5*f(x(i-1),y(i-1))+19*f(x(i),y(i))...
+9*f(x(i+1),s(i+1,j)));
if (s(i+1,j+1)-s(i+1,j))<10^(-6)
y(i+1)=s(i+1,j+1);
break
end
end
end
plot(x,y,'*')
hold on
u=eval(dsolve('Dy=0.09*y-20000','y(0)=200000','x'));
plot(x,u,'r')
end
Comments
Post a Comment