problem in solving. any help will be appreciated.

syms r
kappa = 1; % Specific value for kappa
theta_k = pi/4; % Specific value for theta_k
nu_P = 0.3; % Some value for Poisson's ratio
r_out = 75; % Outer radius
R = 120; % Radius
% Calculate lambda
lambda = (1 + nu_P) / (1 - nu_P);
% Calculate c
c = sqrt(4 - (r_out/R)^2);
% Define the symbolic equations
syms dr(r) dtheta(r);
% Define the equations
dr_eqn = (lambda + 1)*(r*diff(dr, r, r) + diff(dr, r)) + (1/r)*((kappa*r^2*cos(theta_k) - (lambda + 1))*dr) - kappa*r*dtheta*sin(theta_k) == -16*lambda*r^2/(c^4*R^2);
dtheta_eqn = r*diff(dtheta, r, r) + diff(dtheta, r) + (1/r)*((kappa*r^2*cos(theta_k) - 1)*dtheta) + kappa*r*dr*sin(theta_k) == 0;
% Substitute boundary conditions into equations
dr_eqn_subs = subs(dr_eqn, {dr(0), dr(r_out)}, [0, 0]);
dtheta_eqn_subs = subs(dtheta_eqn, {dtheta(0), dtheta(r_out)}, [0, 0]);
% Solve the equations
sol_dr = dsolve(dr_eqn_subs);
sol_theta = dsolve(dtheta_eqn_subs);
% Define parameter values
r_values = linspace(0, r_out, 100); % Corrected range
% Evaluate solutions for specific values of kappa and theta_k
d_r_fun = matlabFunction(subs(sol_dr, {kappa, theta_k}, {kappa, theta_k}));
d_theta_fun = matlabFunction(subs(sol_theta, {kappa, theta_k}, {kappa, theta_k}));
% Plot solutions
figure;
subplot(2, 1, 1);
plot(r_values, d_r_fun(r_values));
xlabel('r');
ylabel('d_r(r)');
title(['d_r(r) vs r for \kappa = ', num2str(kappa), ', \theta_k = ', num2str(theta_k)]);
grid on;
subplot(2, 1, 2);
plot(r_values, d_theta_fun(r_values));
xlabel('r');
ylabel('d_{\theta}(r)');
title(['d_{\theta}(r) vs r for \kappa = ', num2str(kappa), ', \theta_k = ', num2str(theta_k)]);
grid on;

want to numberically solve the eqns and plot d_r(r) vs r and d_theta(r) vs r.

[Sam Chak]

Could you please provide the differential equations? This will enable us to cross-reference them and verify if they are accurately represented in the code.

[tuhin]

these two are the coupled differential eqns. which I want to solve numerically to get plots of d_r(r) vs r and d_theta(r ) vs r. Here one can use kappa and theta_k as variables. Please check and let me know if you need any more details

[tuhin]

I have the analytic solutions of this. Please see the attached. I need to plot d_r(r) vs r and d_theta(r ) vs r. the solutions are:

[Sam Chak]

It is a boundary value problem. Try using the bvp4c() solver. Formulate the problem as shown in the example given in the link.

Answer 1

[SAI SRUJAN]

Hi Tuhin,

I understand that you are facing an issue issue solving the equations related to the boundary value problem.

We can use 'bvp4c' or 'bvp5c' to solve boundary value problems (BVPs) for ordinary differential equations in MATLAB. These functions are particularly suited for problems where you know the conditions at both ends of the interval. Convert the second-order ODEs into a system of first-order ODEs, as 'bvp4c' and 'bvp5c' require the equations to be in first-order form.

Please go through the following code sample to proceed further,

function solveBVP()
    kappa = 1;
    theta_k = pi/4;
    nu_P = 0.3;
    r_out = 75;
    R = 120;
    lambda = (1 + nu_P) / (1 - nu_P);
    c = sqrt(4 - (r_out/R)^2);

    sol = bvp4c(@odefun, @bcfun, solinit);

    r = linspace(0, r_out, 100);
    y = deval(sol, r);

    figure;
    subplot(2, 1, 1);
    plot(r, y(1, :));
    xlabel('r');
    ylabel('d_r(r)');
    title('d_r(r) vs r');
    grid on;

    subplot(2, 1, 2);
    plot(r, y(2, :));
    xlabel('r');
    ylabel('d_{\theta}(r)');
    title('d_{\theta}(r) vs r');
    grid on;
end

function dydr = odefun(r, y)
    dydr = [...]; % System of first-order ODEs derived from your second-order equations
end

function res = bcfun(ya, yb)
    % Boundary conditions at r=0 and r=r_out
    % ya represents the solution at the start of the interval
    % yb represents the solution at the end of the interval
    res = [...];
end

% Initial guess
solinit = bvpinit(linspace(0, r_out, 10), [0 0 0 0]);

We need to fill in the 'odefun' with the system of equations converted to first-order form, and 'bcfun' with the specific boundary conditions. We need to adjust the functions 'odefun','bcfun' and 'solinit' as per the nature of the problem to get the desired output.

For a comprehensive understanding of the 'bvp4c' function in MATLAB, please refer to the following documentation.

https://www.mathworks.com/help/matlab/ref/bvp4c.html

I hope this helps!