Scribd
Upload
3 views2 pages

MATLAB BVP and RK4 Methods Explained

The document contains MATLAB code for solving a boundary value problem (BVP) and implementing numerical methods such as Runge-Kutta and shooting methods. It defines functions for calculating the system's behavior under specific parameters and plots the results for comparison. The code includes nested functions for the numerical integration and defines the system dynamics based on given physical parameters.

Uploaded by

nehar shubhescha
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as TXT, PDF, TXT or read online on Scribd
3 views2 pages

MATLAB BVP and RK4 Methods Explained

The document contains MATLAB code for solving a boundary value problem (BVP) and implementing numerical methods such as Runge-Kutta and shooting methods. It defines functions for calculating the system's behavior under specific parameters and plots the results for comparison. The code includes nested functions for the numerical integration and defines the system dynamics based on given physical parameters.

Uploaded by

nehar shubhescha
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as TXT, PDF, TXT or read online on Scribd

function termfinal()

w=15e3;I=0.0003; L=3; E=200e9;


b=0;
f=@(x,y) [y(2); ((w*L*x/2*E*I)-(w*x.^2/2*E*I))];
boundary_cond=@(ya, yb) [ya(1)-b; yb(1)-b];
sol_initial=bvpinit(linspace(0,3,100),[b b]);
sol=bvp4c(f, boundary_cond,sol_initial);
x = linspace(0,3);
y = deval(sol,x);
plot(x,y(1,:),'ko')
grid on
format longG
y_exact= (((2*w*L*x.^3)-(w*x.^4)-(w*L^3*x))/24*E*I);
hold on
plot(x,y_exact,'r')
hold off
%multi_segment_RK4
w=15e3; L=3; I=30000e-8; E=200e9;
h=0.6;
xspan=[0 3];
n=(xspan(2)-xspan(1))/h;
C=[1 0; 0 0];
D=[0 0;1 0];
F=[0;0];
x=xspan(1):h:xspan(2);

[x, ya]=MyRK4(@dY1dx, xspan, [1 0], h);


[x, ya1]=MyRK4(@dY1dx, xspan, [0 1], h);
ya=ya'
ya1=ya1'
Yap=[ya(:,end) ya1(:,end)];
[x, z1]=MyRK4(@dzdx, xspan, [0 0], h);
z1=z1';
Zp=z1(:,end);
ya01=([C+D*Yap]^-1)*(F-D*Zp);
for i=1:length(x)
Yap1=ya(:,i);
Yap2=ya1(:,i);
Yap12=[Yap1 Yap2];
Zp1=z1(:,i);
yp1=Yap12*ya01+Zp1;
y1(i)=yp1(1);
end
hold on
plot(x,y1,'m-')
grid on
hold off
%% nested function
function [x,y]=MyRK4(odefun,xspan,y0,h)
x=xspan(1):h:xspan(2);
y(1,:)=y0;
for n=1:length(x)-1
k1=odefun(x(n),y(n,:))';
k2=odefun(x(n)+(h/2),y(n,:)+(k1*h/2))';
k3=odefun(x(n)+(h/2),y(n,:)+(k2*h/2))';
k4=odefun(x(n)+h,y(n,:)+k3*h)';
y(n+1,:)=y(n,:)+(k1+2*k2+2*k3+k4)*h/6;
end
end
function xp=dzdx(x,y)
w=15e3; L=3; I=0.0003; E=200e9;
xp(1)=y(2);
xp(2)=((w*x*(L-x)/2*E*I));
xp=xp';
end
function xp=dY1dx(x,y)
xp(1)=y(2);
xp(2)=0;
xp=xp';
end
%Shooting Method
w=15e3;I=0.0003; L=3; E=200e9;
x0 = 0; xb = 3; h = 0.6; y0 = 0; yb = 0;
%when zo=1
z01 = 1;
[x1,y1]=RK4(@dydx,[x0 xb],[y0 z01],h);
disp(y1(end,1));
%when zo=4
z02 = 4;
[x2,y2]=RK4(@dydx,[x0 xb],[y0 z02],h);
disp(y2(end,1));
z03 = z01 + (z02-z01)*(yb-y1(end,1))/(y2(end,1)-y1(end,1));
[x3,y3]=RK4(@dydx,[x0 xb],[y0 z03],h);
disp(y3(:,1));
hold on
plot(x3,y3(:,1),'*', 'LineWidth',1)
grid on
legend('bvp','Exact','multi segment', 'shooting')
title('Term Assignment P01');
hold off
function [x,y]=RK4(ODEfunc,xspan,T0,h)
x=xspan(1):h:xspan(2);
y(1,:)=T0;
for n=1:length(x)-1
k1=ODEfunc(x(n),y(n,:))';
k2=ODEfunc(x(n)+(h/2),y(n,:)+(h/2)*k1)';
k3=ODEfunc(x(n)+(h/2),y(n,:)+(h/2)*k2)';
k4=ODEfunc(x(n)+h,y(n,:)+h*k3)';
y(n+1,:)=y(n,:)+((1/6)*k1+(1/3)*(k2+k3)+(1/6)*k4)*h;
end
end
function yp=dydx(x,y)
yp(1)=y(2);
yp(2)=((w*x*(L-x)/2*E*I));
yp=yp';
end
end

You might also like