############################################################################### ## Solution of Partial Differential Equation by Oberst-Riquier Algorithm ## ## ## ## Author: Ayan Sengupta, Debasattam Pal ## ## Version: 1.0 ## ############################################################################### reset() # Resets all variables previously defined import timeit as tm import re exp = 10 ''' Here Oberst-Riquier algorithm is used for solving the given set of PDE's. The output of this algorithm is in a power series form. The program first asks for the number of independent variables, after that it asks for the partial differential equation(s). After that the initial conditions must be given as specified. The final solution will be saved in a text file in the current directory named 'output.txt'. To see the solution in the terminal type 'solution'. Once the number of independent variables are defined, variables in d i.e. d1,d2,d3...and variables in x i.e. x1,x2,x3... are created as required. Here di signifies partial derivative with respect to xi. Now, the user needs to give the PDE's in polynomial form in terms of d1,d2,d3... EXAMPLE 1: -Number of independent variables: 2 -Defining d1, d2 -Number of partial differential equation(s): 1 -Enter the f(d) in variables of d -f1= d1-d2^2 -groebner basis generated in 0.049 secs -all monomials generated in 0.051 secs -standard monomial calculated in 0.308 secs -Give initial condition of pde in this form: -Poly(x2)*exp(x2) + -Enter the initial condition: x2*e^(5*x2) -evaluation complete -The solution of the PDE's are saved in file named solution.txt -To see the solution in the terminal, type 'solution' EXAMPLE 2: -Number of independent variables: 3 -Defining d1, d2, d3 -Number of partial differential equation(s): 2 -Enter the f(d) in variables of d -f1= d1-d2*d3^2 -f2= d2+d3*d1^2+3 -groebner basis generated in 0.076 secs -all monomials generated in 0.086 secs -standard monomial calculated in 3.706 secs -Give the initial condition in terms of 1, x3, x3^2, x3^3, x3^4, x3^5, x3^6, x3^7, x3^8, x3^9, x3^10, x2, x2*x3, x2*x3^2, x2*x3^3, x2*x3^4, x2*x3^5, x2*x3^6, x2*x3^7, x2*x3^8, x2*x3^9, x2*x3^10, x2*x3^11, x2*x3^12, x2^2, x2^2*x3, x2^2*x3^2, x2^2*x3^3, x2^2*x3^4, x2^3, x2^3*x3, x2^3*x3^2, x2^3*x3^3, x2^3*x3^4, x2^4, x2^4*x3, x2^4*x3^2, x2^4*x3^3, x2^4*x3^4, x2^5, x2^5*x3, x2^5*x3^2, x2^5*x3^3, x2^5*x3^4, x2^6, x2^6*x3, x2^6*x3^2, x2^6*x3^3, x2^6*x3^4, x2^7, x2^7*x3, x2^7*x3^2, x2^7*x3^3, x2^7*x3^4, x2^8, x2^8*x3, x2^8*x3^2, x2^8*x3^3, x2^8*x3^4, x2^9, x2^9*x3, x2^9*x3^2, x2^9*x3^3, x2^9*x3^4, x2^10, x2^10*x3, x2^10*x3^2, x2^10*x3^3, x2^10*x3^4, x2^11, x2^11*x3, x2^11*x3^2, x2^11*x3^3, x2^11*x3^4, x2^12, x2^12*x3, x2^12*x3^2, x2^12*x3^3, x2^12*x3^4, x2^13, x2^13*x3, x2^13*x3^2, x2^13*x3^3, x2^13*x3^4, x2^14, x2^14*x3, x2^14*x3^2, x2^14*x3^3, x2^14*x3^4, x2^15, x2^15*x3, x2^15*x3^2, x2^15*x3^3, x2^15*x3^4, x2^16, x2^16*x3, x2^16*x3^2, x2^16*x3^3, x2^16*x3^4... -Enter the initial condition: e^5*x2 + (x2+x2^2)*e^x3 -evaluation complete -The solution of the PDE's are saved in file named solution.txt -To see the solution in the terminal, type 'solution' HOW TO RUN THE CODE IN SAGE: First open SAGE from the terminal. Then you can load a file in a Sage session with the command "load". example: if your file is in the current directory where you started sage then: -load("your_file.sage") if your file is in another directory then : -load("~/the_file_location/your_file.sage") ''' def Polynomial_ring(num, precision, mono_order='lex'): ''' This Function defines a polynomial ring with variables d1,d2,d3.... The number of variables defined will be given by the user''' RF = RealField(precision) P = PolynomialRing( RR, ['d%s' %p for p in range(1, num + 1)], order=mono_order) P.inject_variables() return P def Groebner(num): '''This function takes in the polynomials given by the user, generates the ideal and finally computes the reduced groebner basis''' global start_time poly = [0] * num for i in range(num): print 'f' + str(i + 1) + '=', input_temp = str(raw_input()).strip(" ") temp = re.sub('\^','**',input_temp) poly[i] = eval(temp) start_time = tm.default_timer() I = ideal(poly) groebner = I.groebner_basis() return groebner, poly def generate_all_monomial(ring,exp): '''generates all the possible monomials with max degree of each variable defined by the user, d1 d2 with exp3 will generate 15 terms d1,d1^2,d1^3,d*1d2,d1^2*d2...d1^3*d1^3''' var_list = ring.gens_dict().values() num = [exp + 1] * len(var_list) d_v = monomials(var_list, num) d_v.remove(1) return d_v def standard_monomial(all_monomial_list): '''This function generates a set of standard monomials''' T = dict() T_values = [] T_keys = [] for i in range(len(all_monomial_list)): if all_monomial_list[i].reduce(groebner) != 0: T[all_monomial_list[i]] = all_monomial_list[i].reduce(groebner) T_keys = T.keys() T_values = T.values() T_values_coeff = [T_values[i].coefficients()[0] for i in range(len(T_values))] std_monomial_list = [set(T_values[i].monomials()) for i in range(len(T_values))] std_monomial_set = set().union(*std_monomial_list) std_monomial = sorted(list(std_monomial_set)) return T, T_values, T_keys, T_values_coeff, std_monomial def standard_monomial_evaluated_over_initial_condition( initial_condition, std_monomial): std_mon_coeff_from_initial_condition = {} for i in range(len(std_monomial)): derivative_order = list(std_monomial[i].degrees()) diff_parameter = [[x_var[k]] * derivative_order[k] for k in range(len(derivative_order)) if derivative_order[k] != 0] flat_diff_parameter = [ item for sublist in diff_parameter for item in sublist] coeff_temp = initial_condition.derivative(flat_diff_parameter) for p in range(len(x_var)): m = x_var[p] coeff_temp = coeff_temp.subs(m == 0) std_mon_coeff_from_initial_condition[std_monomial[i]] = coeff_temp return std_mon_coeff_from_initial_condition def generate_var_x(num): '''Defines the variables in terms of x1, x2, x3, ... ''' x_var = [var('x%d' % i) for i in range(1, num + 1)] return x_var def generate_factorial(monomial): fac_temp = 1 for p in range(len(monomial.degrees())): fac_temp = fac_temp * factorial(monomial.degrees()[p]) return fac_temp def coeff_temp(f, n): k = len(value[0].degrees()) v = x_var for i in range(k): z = v[i] g = f.derivative(z, value[n].degrees()[i]) f = g.subs(z == 0) return f def calculate_a0(f): for i in range(len(value[0].degrees())): z = x_var[i] f = f.subs(z == 0) return f def generate_solution( std_mon_coeff_from_initial_condition, initial_condition, prec, x_var,approx = 4): coeff_temp = [] for i in range(len(value)): sub_monomial = value[i].monomials() sub_coefficient = value[i].coefficients() sub_monomial_value = [ std_mon_coeff_from_initial_condition[ sub_monomial[k]] for k in range( len(sub_monomial))] total_sub_monomial_sum = sum( [(a * b) for a, b in zip(sub_coefficient, sub_monomial_value)]) coeff_temp.append( float(total_sub_monomial_sum) / generate_factorial( key[i])) #solution_temp = sum([(a * b) for a, b in zip(coeff_temp, key)]) #solution_temp = solution_temp + calculate_a0(initial_condition) solution_temp_approx = sum([(round(a, approx) * b) for a, b in zip(coeff_temp, key)]) solution_temp_approx = solution_temp_approx + \ calculate_a0(initial_condition) solution_approx = solution_temp_approx.subs( {poly_ring.gens()[i]: x_var[i] for i in range(0, num_var)}) solution_temp = [(a * b) for a, b in zip(coeff_temp, key)] solution_temp.append(calculate_a0(initial_condition)) solution_sum = sum(solution_temp) solution = solution_sum.subs( {poly_ring.gens()[i]: x_var[i] for i in range(0, num_var)}) return solution, solution_approx def pde_solution_check(poly, solution): return nil def print_initial_condition(std_monomial, exp): if len(x_var) == 2: value_for_infinity = exp - 2 std_monomial_deg = [] for i in range(len(std_monomial)): std_monomial_deg.append(std_monomial[i].degrees()) temp = std_monomial_deg for i in range(len(std_monomial_deg[0])): for j in range(len(std_monomial_deg)): if std_monomial_deg[j][i] == value_for_infinity: a = std_monomial_deg[j][len(std_monomial_deg[0]) - 1 - i] temp = [value for value in temp if value[ len(std_monomial_deg[0]) - 1 - i] != a] if a == 0: print 'Poly({0})*exp({1}) +'.format(x_var[i], x_var[i]), elif a == 1: print 'Poly({0})*exp({1}) +'.format(x_var[len(std_monomial_deg[0]) - 1 - i], x_var[i]), else: print 'Poly({0}^{1})*exp({2}) +'.format(x_var[len(std_monomial_deg[0]) - 1 - i], std_monomial_deg[j][len(std_monomial_deg[0]) - 1 - i], x_var[i]), #print temp for i in range(len(temp)): if temp[i][0] == 0 and temp[i][1] == 0: print 'Constant + ', elif temp[i][0] == 0: print 'Poly({0}^{1}) +'.format(x_var[1], temp[i][1]), elif temp[i][1] == 0: print 'Poly({0}^{1}) +'.format(x_var[0], temp[i][0]), else: print 'Poly({0}^{1}*{2}^{3}) +'.format(x_var[0], temp[i][0], x_var[1], temp[i][1]), ############################################################# #################### MAIN PROGRAM ########################### ############################################################# num_var = input("\nNumber of independent variables: ") #o = raw_input("monomial ordering(lex, deglex, degrevlex ...): ") poly_ring = Polynomial_ring(num_var, 20) p = input("Number of partial differential equation(s): ") print ("\nEnter the f(d) in variables of d") groebner, poly = Groebner(p) time_groebner = tm.default_timer() print('\ngroebner basis generated in {0} secs'.format( round(time_groebner - start_time, 3))) all_monomial_list = generate_all_monomial(poly_ring,exp) time_allmonomials = tm.default_timer() print('all monomials generated in {0} secs'.format( round(time_allmonomials - start_time, 3))) x_var = generate_var_x(num_var) dic, value, key, coeff, std_monomial = standard_monomial(all_monomial_list) if len(dic) == 0: print ('w = 0') time_std_monomial = tm.default_timer() print('standard monomial calculated in {0} secs\n'.format( round(time_std_monomial - start_time, 3))) #print ('standard monomial set is {0}' .format(std_monomial)) init_cond_temp = str(std_monomial).replace('[','').replace(']','')+'...' init_cond = re.sub('d','x',init_cond_temp) ############################################################# ##################### Initial condition ##################### ############################################################# if len(x_var) == 2: print "\nGive initial condition of pde in this form:\n", else: print("\nGive the initial condition in terms of: {0}".format(init_cond)) print_initial_condition(std_monomial, exp) i = raw_input("\n\nEnter the initial condition: ") f(x) = i print '\nevaluation complete' monomial_coeff = standard_monomial_evaluated_over_initial_condition( f(x), std_monomial) solution, solution_approx = generate_solution(monomial_coeff, f(x), 20, x_var) ############################################################## ####################### OUTPUT ############################# ############################################################## print("\nThe solution of the PDE's are saved in file named solution.txt") print("To see the solution in the terminal, type 'solution'") #print(solution) filepath = 'pde_output.txt' f = open(filepath, 'w') f.write(str(solution)) f.close()