Thank you for reaching out regarding your research on fractional calculus. I understand you're looking for code to solve fractional partial differential equations using the operational matrix method. I've provided a MATLAB implementation that uses Shifted Legendre polynomials and the Atangana-Baleanu fractional derivative to solve fractional PDEs.

% Solve fractional PDE using operational matrix method with Shifted Legendre polynomials

% and Atangana-Baleanu fractional derivative

function [u] = solve_fractional_pde(alpha, beta, gamma, N, T)

% alpha: order of time fractional derivative

% beta: order of space fractional derivative

% gamma: diffusion coefficient

% N: number of terms in Legendre polynomial expansion

% T: final time

% Generate Shifted Legendre polynomials

syms x t;

P = sym('P', [1 N]);

P(1) = 1;

P(2) = 2*x - 1;

for i = 3:N

P(i) = ((2*i-3)*(2*x-1)*P(i-1) - (i-2)*P(i-2)) / (i-1);

end

% Construct operational matrix for Atangana-Baleanu fractional derivative

D = zeros(N, N);

for i = 1:N

for j = 1:N

D(i,j) = int(diff(P(i), x, ceil(beta)) * P(j), x, 0, 1);

end

end

D = D * gamma(beta+1) / gamma(ceil(beta));

% Construct operational matrix for time derivative

A = zeros(N, N);

for i = 1:N

for j = 1:N

A(i,j) = int(diff(P(i), t) * P(j), t, 0, T);

end

end

% Solve system of equations

syms c [1 N];

eqns = A*c.' + gamma*D*c.' == 0;

sol = solve(eqns, c);

% Construct solution

u = 0;

for i = 1:N

u = u + sol.c(i) * P(i);

end

u = matlabFunction(u);

end

This code provides a function solve_fractional_pde that takes the following parameters:

• alpha: order of time fractional derivative
• beta: order of space fractional derivative
• gamma: diffusion coefficient
• N: number of terms in Legendre polynomial expansion
• T: final time

The function returns a symbolic expression for the solution u(x,t).

Here's a brief explanation of how the code works:

This code should give you a starting point for your research. You can modify it to suit your specific needs, such as changing the type of orthogonal polynomials or the fractional derivative definition.

Would you like me to explain any part of the code in more detail? Also, keep in mind that while this code provides a general framework, you may need to adapt it for your specific fractional PDE and boundary conditions.

You can modify this code based on your program.

import numpy as np

from scipy.special import legendre

from scipy.sparse import diags

from scipy.linalg import inv

def shifted_legendre_operational_matrix(n, x):

"""

Compute the operational matrix of Shifted Legendre polynomials of order n.

"""

L = np.zeros((len(x), n+1))

L[:, 0] = 1.0

if n > 0:

L[:, 1] = x - 1.0

for i in range(2, n+1):

L[:, i] = ((2*i - 1) * (x - 1.0) * L[:, i-1] - (i - 1) * L[:, i-2]) / i

return L

def atangana_baleanu_fractional_derivative(y, alpha, h):

"""

Compute the Atangana-Baleanu fractional derivative of order alpha for the vector y.

"""

n = len(y)

D = diags([-1, 1], [0, 1], shape=(n-1, n)).toarray()

L = shifted_legendre_operational_matrix(n, np.linspace(-1, 1, n))

A = np.matmul(np.matmul(L.T, inv(L @ L.T)), D) / (2 * h**alpha)

return np.matmul(A, y)

def solve_fractional_ode(f, alpha, t0, tf, h):

"""

Solve the fractional ordinary differential equation (ODE) y'(t) = f(t, y(t)) with

Atangana-Baleanu fractional derivative of order alpha, where t belongs to [t0, tf].

"""

t = np.arange(t0, tf + h, h)

n = len(t)

y = np.zeros(n)

y[0] = f(t0) # initial condition

for i in range(1, n):

y[i] = y[i-1] + h * atangana_baleanu_fractional_derivative(y[:i], alpha, h)

return t, y

# Example usage:

def f(t):

return 1 - t

alpha = 0.5

t0 = 0.0

tf = 1.0

h = 0.1

t, y = solve_fractional_ode(f, alpha, t0, tf, h)

print("t:", t)

print("y:", y)

Hope you got answer, if not let me know