Scribd
Upload
127 views

Python

This document contains 8 questions involving Python code to solve various mathematical problems using libraries like NumPy, SciPy, and Matplotlib. Question 1 animates optimization of a 2D function using conjugate gradient and BFGS methods. Question 2 optimizes the positions of 13 atoms to minimize potential energy. Question 3 performs linear interpolation of sin(x) at 4 points. Question 4 fits a 3rd degree polynomial to noisy sin(x/2) data. Question 5 numerically integrates 1/(1+x^2) from 0 to infinity. Question 6 performs double integration of an exponential function over a region. Question 7 solves an ODE initial value problem. Question 8 solves a system of 2 coupled ODEs.

Uploaded by

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

Python

This document contains 8 questions involving Python code to solve various mathematical problems using libraries like NumPy, SciPy, and Matplotlib. Question 1 animates optimization of a 2D function using conjugate gradient and BFGS methods. Question 2 optimizes the positions of 13 atoms to minimize potential energy. Question 3 performs linear interpolation of sin(x) at 4 points. Question 4 fits a 3rd degree polynomial to noisy sin(x/2) data. Question 5 numerically integrates 1/(1+x^2) from 0 to infinity. Question 6 performs double integration of an exponential function over a region. Question 7 solves an ODE initial value problem. Question 8 solves a system of 2 coupled ODEs.

Uploaded by

Kalim Ullah Marwat
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
You are on page 1/ 10

Python Assignment

Name : Kalim ullah Chinese name: wu kailin


student ID: LX 20230512103

Laoshi: 赛欧拉

College of Science
Department mathematics
Question No 1
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from matplotlib.animation import FuncAnimation

# Define the function


def f(x):
return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

# Define the gradient of the function


def df(x):
dfdx0 = -2 * (1 - x[0]) - 400 * x[0] * (x[1] - x[0]**2)
dfdx1 = 200 * (x[1] - x[0]**2)
return np.array([dfdx0, dfdx1])

# Initial guess
initial_guess = np.array([0, 0])

# Optimize using Conjugate Gradient (CG) method


result_cg = minimize(f, initial_guess, method='CG', jac=df, options={'disp': True})

# Optimize using Broyden–Fletcher–Goldfarb–Shanno (BFGS) method


result_bfgs = minimize(f, initial_guess, method='BFGS', jac=df, options={'disp': True})

# Plot the animated figures


fig, ax = plt.subplots(1, 2, figsize=(12, 5))
ax[0].set_title('Conjugate Gradient (CG) Method')
ax[1].set_title('BFGS Method')

# Function contour
x_vals = np.linspace(-2, 2, 400)
y_vals = np.linspace(-1, 3, 400)
X, Y = np.meshgrid(x_vals, y_vals)
Z = f([X, Y])
contour = ax[0].contour(X, Y, Z, levels=np.logspace(-1, 3, 10), cmap='viridis')
contour = ax[1].contour(X, Y, Z, levels=np.logspace(-1, 3, 10), cmap='viridis')

# Plot the optimization path


def plot_optimization_path(ax, result, method_name):
x_opt = result.x
ax.plot(x_opt[0], x_opt[1], 'ro', label='Optimal Point', markersize=8)
ax.legend()

# Animation update function


def update(frame):
ax[0].cla()
ax[1].cla()
ax[0].set_title('Conjugate Gradient (CG) Method')
ax[1].set_title('BFGS Method')
ax[0].contour(X, Y, Z, levels=np.logspace(-1, 3, 10), cmap='viridis')
ax[1].contour(X, Y, Z, levels=np.logspace(-1, 3, 10), cmap='viridis')
plot_optimization_path(ax[0], result_cg, 'Conjugate Gradient')
plot_optimization_path(ax[1], result_bfgs, 'BFGS')

# Create the animation


animation = FuncAnimation(fig, update, frames=range(100), repeat=False)
plt.show()
Question No 2 There 13 atoms, for each pair of atoms, there is a potential r ^ a * (- 6) - r ^ a * (- 12) The
total potential of these 13 atoms are the potential sum of each pair. Find the global minima potential
of these 13atoms.

import numpy as np
from scipy.optimize import minimize

# Number of atoms
n_atoms = 13

# Function to calculate the potential energy for a pair of atoms


def potential_pair(r, s):
return r**s * (-6) - r**s * (-12) # Corrected the subtraction

# Function to calculate the total potential energy for all pairs of atoms
def total_potential(positions):
energy = 0
for i in range(n_atoms): # Corrected the variable name from I to i
for j in range(i + 1, n_atoms): # Corrected the variable name from J to j
r = np.linalg.norm(positions[i] - positions[j])
energy += potential_pair(r, 2) # You can adjust the value of ‘s’ as needed

return energy

# Initial random positions for optimization


initial_positions = np.random.rand(n_atoms * 3) # Flatten to one-dimensional array

# Define the optimization problem


result = minimize(total_potential, initial_positions, method='L-BFGS-B') # Corrected the
quotation marks
# Extract the optimized positions and energy
optimized_positions = result.x.reshape((n_atoms, 3))
global_min_potential = result.fun # Corrected the variable name from Global_min_potential to
global_min_potential
print("Optimized Positions:")
print(optimized_positions)
print("Global Minimum Potential:", global_min_potential

Question No 3 interplate f(x) = sin(x) on x = 0 ,pi /2, pi, 3* pi / 2

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

# Define the function f(x) = sin(x)


def f(x):
return np.sin(x)

# Define the interpolation points


x_points = np.array([0, np.pi/2, np.pi, 3*np.pi/2])

# Compute the function values at the interpolation points


y_points = f(x_points)

# Create the interpolation function


interpolated_function = interp1d(x_points, y_points, kind='linear')

# Define a range of x values for plotting the original function and the interpolation
x_values = np.linspace(0, 3*np.pi/2, 1000)
y_original = f(x_values)
y_interpolated = interpolated_function(x_values)

# Plot the original function and the interpolation


plt.plot(x_values, y_original, label='Original function: $sin(x)$')
plt.scatter(x_points, y_points, color='red', marker='o', label='Interpolation Points')
plt.plot(x_values, y_interpolated, linestyle='dashed', label='Linear Interpolation')
plt.legend()
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Interpolation of $sin(x)$')
plt.show()
Question 4: Fit sin(x/2) on interval [- 1, 1] with polynomial of degree 3

import numpy as np
import matplotlib.pyplot as plt

# Define the true function


def true_function(x):
return np.sin(x/2)

# Generate data points for fitting


np.random.seed(42) # For reproducibility
x_data = np.linspace(-1, 1, 100)
y_data = true_function(x_data) + np.random.normal(0, 0.1, len(x_data))

# Fit a polynomial of degree 3


degree = 3
coefficients = np.polyfit(x_data, y_data, degree)

# Generate x values for plotting the true function and the fitted polynomial
x_values = np.linspace(-1, 1, 1000)
y_true = true_function(x_values)
y_fit = np.polyval(coefficients, x_values)

# Plot the true function, noisy data points, and the fitted polynomial
plt.plot(x_values, y_true, label='True Function')
plt.scatter(x_data, y_data, color='red', label='Noisy Data Points')
plt.plot(x_values, y_fit, label=f'Polynomial Fit (Degree {degree})', linestyle='dashed')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()
Question 5 integration

from scipy import integrate

# Define the integrand function


def integrand(x):
return 1 / (1 + x**2)

# Perform the numerical integration from 0 to infinity


result, error = integrate.quad(integrand, 0, float('inf'))

# Print the result


print("The result of the integration from 0 to infinity is:", result)
Question No 6

import numpy as np
from scipy.integrate import dblquad

# Define the integrand


def integrand(x, y):
return np.exp(-y**2)

# Define the limits of integration


def x_lower_limit(y):
return 0

def x_upper_limit(y):
return y

# Perform the double integration


result, error = dblquad(integrand, 0, 1, x_lower_limit, x_upper_limit)

print("Result of the double integral:", result)


Question 7 Solve (x ^ 2 – 1) * y’ + 2x * y ^ 2 = 0’ with initial value y(2) = 1
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Define the ODE as a function


def model(y, x):
return (1 - x**2) * y / (2 * x)

# Initial condition
y0 = 1
# Time points for integration
x = np.linspace(2, 10, 100) # Adjust the time range as needed

# Solve the ODE


solution = odeint(model, y0, x)
# Extract the solution for y
y_values = solution[:, 0]

# Plot the result


plt.plot(x, y_values, label='Solution')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of the ODE: (1 – x^2) * y / (2 * x)')
plt.legend()
plt.show()
Question No 8

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Define the system of differential equations


def model(u, t):
x, y = u
dxdt = 2*x + y
dydt = 3*x + 2*y
return [dxdt, dydt]

# Initial conditions
u0 = [2, 0]

# Time points for integration


t = np.linspace(0, 5, 100) # Adjust the time range as needed

# Solve the system of differential equations


solution = odeint(model, u0, t)

# Extract the solutions for x and y


x_values = solution[:, 0]
y_values = solution[:, 1]

# Plot the results


plt.plot(t, x_values, label='x(t)')
plt.plot(t, y_values, label='y(t)')
plt.xlabel('time')
plt.ylabel('Values')
plt.legend()
plt.show()

You might also like