MATLAB实现二重积分计算

MATLAB实现二重积分计算

太原理工大学中外合作办学-数值分析作业题

1. Question 1: (50 marks)

Evaluate the following double integral. Note that the upper y boundary is a function of x

(a)Manually solve this question by using one Simpson 1/3 rule in the x-direction and two Simpson 1/3 rules in the y-direction. (15 marks)
(It is requested to show the detailed calculation procedure in your report, including the intervals, meshing, function values at each node, and integral for x and y, etc)
(b)Develop a MATLAB M-file based on the Simpson 1/3 rules for evaluation of the double integral of the above function with any number of Simpson 1/3 rules in the x- and y-directions. Note the number of Simpson 1/3 may be different in the x- and y-directions

Then, use the developed code to calculate the above double integral by different numbers of Simpson 1/3 rules until an accurate result is achieved. (25 marks)

(The M-file must be well commented and included in the report, and the developed code must be submitted for checking. You need to try different numbers of Simpson 1/3 rules in x- and y-directions until an accurate result is achieved. And justify why the final result is acceptable).

©Choose a proper MATLAB built-in function and develop an M-file script to calculate the above double integral.
(10 marks)
(The M-file must be well commented and included in the report, and the developed code must be submitted for checking. You are requested to compare the results by built-in function with the data obtained by the developed code. )

• 问题1解答

（a）Normal dual integration




Simpson 1/3rule：

（b） 解答

The code implementation using matlab is as follows.

% Define the integrand
f = @(x, y) x^2 - 2*y^2 + x*y^3;
 
% Define the integration region
x_start = 0;
x_end = 2;
y_start = -1;
y_end = @(x) 0.5 + 0.5*x;
 
% Define the number of intervals for the Simpson 1/3 rule in the x and y directions
num_intervals_x = 1;
num_intervals_y = 2;
 
% Calculate the step sizes in the x and y directions
h_x = (x_end - x_start) / num_intervals_x;
 
% Initialize the integration result
result = 0;
 
% Perform double integration using Simpson's 1/3 rule
for i = 1:num_intervals_x
    % Calculate the start and end values of the current x interval
    x1 = x_start + (i-1) * h_x;
    x2 = x_start + i * h_x;
    
    % Calculate the number of intervals for the Simpson 1/3 rule in the y direction for the current x interval
    num_intervals_y = 10 + i;
    
    % Calculate the step size in the y direction
    h_y = (y_end(x2) - y_start) / num_intervals_y;
    
    % Use Simpson's 1/3 rule to calculate the double integral for the current x interval
    for j = 1:num_intervals_y
        % Calculate the start and end values of the current y interval
        y1 = y_start + (j-1) * h_y;
        y2 = y_start + j * h_y;
        
        % Calculate the area of the current region
        area = (x2 - x1) * (y2 - y1);
        
        % Calculate the function values at the boundary points
        f_x1y1 = f(x1, y1);
        f_x2y1 = f(x2, y1);
        f_x1y2 = f(x1, y2);
        f_x2y2 = f(x2, y2);
        f_xmidy1 = f(x1 + (x2 - x1) / 2, y1);
        f_xmidy2 = f(x1 + (x2 - x1) / 2, y2);
        
        % Use Simpson's 1/3 rule to calculate the integral contribution of the current region
        integral_contribution = (area/9) * (f_x1y1 + f_x2y1 + f_x1y2 + f_x2y2 + 4*f_xmidy1 + 4*f_xmidy2);
        
        % Accumulate the integral contribution of the current region to the result
        result = result + integral_contribution;
    end
end
 
% Display the result
disp(result);

By choosing different num_intervals_x and num_intervals_y comparisons, it is found that when num_intervals_x=1 and num_intervals_y=2, the result is 3.7477, which is closest to 3.6396.
The above code runs as shown below.

（c）解答
The built-in function integral2 is chosen to implement the dual integration solution, the code is as follows.

% Define the function to be integrated
f = @(x, y) x.^2 - 2*y.^2 + x.*y.^3;
 
% Define the limits for x
x_lower = 0;
x_upper = 2;
 
% Define the lower limit for y and the upper limit as a function of x
y_lower = -1;  % Lower limit for y
y_upper = @(x) 0.5 + 0.5 * x;  % Upper limit for y as a function of x
 
% Use integral2 to calculate the double integral
result = integral2(f, x_lower, x_upper, y_lower, y_upper);
 
disp(['The double integral result using integral2 with variable upper limit is: ', num2str(result)]);

The above code runs with a result of 3.6396, as shown below.

2. Question 2: (50 marks)
The transverse deflection of a cable of length , fixed at both ends, is given as a solution to

where,
T = tension in cable
R = flexural stiffness
q = distributed transverse load

(Three shooting need to be conducted and plot the results; The M-file must be well commented and included in the report, and the developed code must be submitted for checking).

• 问题2解答
  （a）
  
  
  
  （b）

Step1.assume



Step2.create a new variable

Step3.solve by 2nd RK Heun method

Predict u1 and u2:

Hence,x=0.3,u1=0.0009,u2=0.153

©

% Parameters
L = 1.2;
T = 9000;
q = 10;
R = 2e5;
step_size = 0.3;
 
% Initial conditions
x = 0;
u = 0;
v = 0.003;
 
% Arrays to store values for plotting
x_values = [];
u_values = [];
 
% Iteration
while x <= 0.6
    % RK Heun method
    k1u = step_size * v;
    k1v = step_size * (T*u/R + q*x*(x-L)/(2*R));
    
    k2u = step_size * (v + k1v);
    k2v = step_size * (T*(u + k1u)/R + q*(x + step_size)*(x+step_size-L)/(2*R));
    
    % Update position using linear interpolation
    u_temp = u + 0.5 * (k1u + k2u);
    
    % Display result before refining
    fprintf('Before refining - x = %.2f, u = %.4f\n', x + step_size, u_temp);
    
    % Refine using linear interpolation
    v_temp = v + 0.5 * (k1v + k2v);
    k1_refine = step_size * v_temp;
    k2_refine = step_size * (T*u_temp/R + q*(x + step_size)*(x+step_size-L)/(2*R));
    
    % Linear interpolation to refine position
    u = u_temp + 0.5 * (k1_refine + k2_refine);
    
    % Display refined result
    fprintf('After refining - x = %.2f, u = %.4f\n', x + step_size, u);
    
 
    
    % Update velocity for the next iteration
    v = v_temp;
    
    % Update position for the next iteration
    x = x + step_size;
    
    % Store values for plotting
    x_values = [x_values, x];
    u_values = [u_values, u];
end
 
% Plot transverse deflection versus cable length
figure;
plot(x_values, u_values, '-o');
xlabel('Cable Length (L)');
ylabel('Transverse Deflection (u)');
title('Transverse Deflection vs. Cable Length (Shooting Method)');
grid on;

The above code runs as follows.

The relationship between u and L is shown below.

When du/dx=0.004 the implementation code is as follow.

% Parameters
L = 1.2;
T = 9000;
q = 10;
R = 2e5;
step_size = 0.3;
 
% Initial conditions
x = 0;
u = 0;
v = 0.004;
 
% Arrays to store values for plotting
x_values = [];
u_values = [];
 
% Iteration
while x <= 0.6
    % RK Heun method
    k1u = step_size * v;
    k1v = step_size * (T*u/R + q*x*(x-L)/(2*R));
    
    k2u = step_size * (v + k1v);
    k2v = step_size * (T*(u + k1u)/R + q*(x + step_size)*(x+step_size-L)/(2*R));
    
    % Update position using linear interpolation
    u_temp = u + 0.5 * (k1u + k2u);
    
    % Display result before refining
    fprintf('Before refining - x = %.2f, u = %.4f\n', x + step_size, u_temp);
    
    % Refine using linear interpolation
    v_temp = v + 0.5 * (k1v + k2v);
    k1_refine = step_size * v_temp;
    k2_refine = step_size * (T*u_temp/R + q*(x + step_size)*(x+step_size-L)/(2*R));
    
    % Linear interpolation to refine position
    u = u_temp + 0.5 * (k1_refine + k2_refine);
    
    % Display refined result
    fprintf('After refining - x = %.2f, u = %.4f\n', x + step_size, u);
    
 
    
    % Update velocity for the next iteration
    v = v_temp;
    
    % Update position for the next iteration
    x = x + step_size;
    
    % Store values for plotting
    x_values = [x_values, x];
    u_values = [u_values, u];
end
 
% Plot transverse deflection versus cable length
figure;
plot(x_values, u_values, '-o');
xlabel('Cable Length (L)');
ylabel('Transverse Deflection (u)');
title('Transverse Deflection vs. Cable Length (Shooting Method)');
grid on;

The above code runs as follows.

The relationship between u and L is shown below.

（d）

When du/dx=0.003 , this is achieved using the built-in function ode45 with the following code.

% Parameters
L = 1.2;
T = 9000;
q = 10;
R = 2e5;
 
% Define the ODE as a function
ode = @(x, y) [y(2); T*y(1)/R + q*x*(x-L)/(2*R)];
 
% Initial conditions
x0 = 0;
y0 = [0; 0.003]; % [u; v]
 
% Define the integration range
x_span = [0:0.3:0.6]; % Ensure step size is 0.3
 
% Solve the ODE using ode45
[x, y] = ode45(ode, x_span, y0);
 
% Extract the solution
u_solution = y(:, 1);
 
% Display the results
for i = 1:length(x)
    fprintf('x = %.2f, u = %.4f\n', x(i), u_solution(i));
end

The above code runs as follows.

When du/dx=0.0004 , this is achieved using the built-in function ode45 with the following code.

% Parameters
L = 1.2;
T = 9000;
q = 10;
R = 2e5;
 
% Define the ODE as a function
ode = @(x, y) [y(2); T*y(1)/R + q*x*(x-L)/(2*R)];
 
% Initial conditions
x0 = 0;
y0 = [0; 0.004]; % [u; v]
 
% Define the integration range
x_span = [0:0.3:0.6]; % Ensure step size is 0.3
 
% Solve the ODE using ode45
[x, y] = ode45(ode, x_span, y0);
 
% Extract the solution
u_solution = y(:, 1);
 
% Display the results
for i = 1:length(x)
    fprintf('x = %.2f, u = %.4f\n', x(i), u_solution(i));
end

The above code runs as follows.