matlab code for 2nd order differential equation with boundary conditions
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
(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 Kommentare
Antworten (1)
Torsten
am 7 Nov. 2023
Bearbeitet: Torsten
am 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 Kommentare
Siehe auch
Kategorien
Mehr zu Calculus finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!