matlab code for 2nd order differential equation with boundary conditions
6 views (last 30 days)
Show older comments
(d^2 u(x))/〖dx〗^2 =(γu(x))/(1+αu(x))
at x=0,u=1
at x=1,du/dx=0
α=0.001,γ=1,2.5,5,10,100 plot u verses x at differnet γ values
% Parameters
alpha = 0.001;
gamma = [1, 2.5, 5, 10, 100];
% Define the differential equation
du2_dx2 = @(x, u) (gamma * u) / (1 + alpha * u);
% Define the range of x values
x_span = [0, 1]; % from x=0 to x=1
% Initial conditions at x=0
x0 = 0;
u0 = 1;
% Boundary conditions at x=1 (du/dx = 0)
boundary_condition = 0;
% Define function to compute residuals for boundary condition
boundary_residuals = @(u_end) (u_end(1) - boundary_condition);
% Solve the ODE with boundary condition
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[x, u] = ode45(@(x, u) [u(2); du2_dx2(x, u(1))], x_span, [u0, 0], options);
% Adjust u(1) to satisfy the boundary condition du/dx(1) = 0
u_end = u(end, :);
u_end(2) = u_end(2) + boundary_residuals(u_end);
[t, u] = ode45(@(x, u) [u(2); du2_dx2(x, u(1))], x_span, u_end, options);
% Plot u versus x
plot(t, u(:, 1), 'b-', 'LineWidth', 0.2);
xlabel('x');
ylabel('u(x)');
title('Plot of u(x) versus x');
axis([0 1 0 1]); % Set axis limits for x and u
0 Comments
Answers (1)
Torsten
on 7 Nov 2023
Edited: Torsten
on 7 Nov 2023
You have a boundary value problem, not an initial value problem. ode45 is not suited in this case. Use bvp4c or bvp5c instead.
% Parameters
alpha = 0.001;
gamma = [1, 2.5, 5, 10, 100];
hold on
for i = 1:numel(gamma)
fcn = @(x,y)[y(2);gamma(i)*y(1)/(1+alpha*y(1))];
bc = @(ya,yb)[ya(1)-1;yb(2)];
guess = @(x)[1;0];
xmesh = linspace(0,1,20);
solinit = bvpinit(xmesh, guess);
sol = bvp4c(fcn, bc, solinit);
plot(sol.x,sol.y(1,:))
end
hold off
grid on
0 Comments
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!