In Maple 2021, it is now possible to use object:-method(arg) notation. This makes is
easier to use OOP in maple. To do this, use _self as follows
restart; module person() option object; local name::string:=""; local age::integer:=0; export get_name::static:=proc(_self,$) return _self:-name; end proc; export set_name::static:=proc(_self,name::string,$) _self:-name:=name; end proc; export get_age::static:=proc(_self,$) return _self:-age; end proc; export set_age::static:=proc(_self,age::integer,$) _self:-age:=age; end proc; end module;
And now make an object and use it as follows
o:=Object(person) o := Object<<1846759887808>> o:-get_name(); "" o:-get_age(); 0 o:-set_name("joe doe"); o:-get_name(); "joe doe"
Add ModuleCopy proc in the class. This will automatically be called to initialize the
object.
Here is an example
restart; module ODE() option object; local ode:=NULL; local y::symbol; local x::symbol; local sol; export ModuleCopy::static := proc( _self::ODE, proto::ODE, ode, func, $ ) print("Initilizing object with with args: ", [args]); _self:-ode:= ode; _self:-y:=op(0,func); _self:-x:=op(1,func); end proc; export dsolve::static:=proc(_self,$) _self:-sol := :-dsolve(ode,y(x)); end proc; export get_sol::static:=proc(_self,$) return sol; end proc; end module;
And now make an object and use it as follows
o:=Object(ODE, diff(y(x),x)+y(x)=sin(x), y(x)); o:-dsolve(): o:-get_sol(); #y(x) = -1/2*cos(x) + 1/2*sin(x) + exp(-x)*_C1
So a constructor just makes it easier to initialize the object without having to make a
number of set() calls to initialize each memeber data.
This is done using overload with different ModuleCopy proc in the class.
Here is an example. Lets make a constructor that takes an ode and initial conditions, and one that only takes an ode with no initial conditions.
restart; module ODE() option object; local ode:=NULL; local y::symbol; local x::symbol; local ic:=NULL; local sol; export ModuleCopy::static:= overload( [ proc( _self::ODE, proto::ODE, ode, func, $ ) option overload; _self:-ode:= ode; _self:-y:=op(0,func); _self:-x:=op(1,func); end proc, proc( _self::ODE, proto::ODE, ode, func, ic, $ ) option overload; _self:-ode:= ode; _self:-y:=op(0,func); _self:-x:=op(1,func); _self:-ic :=ic; end proc ] ); export dsolve::static:=proc(_self,$) if evalb(ic=NULL) then sol := :-dsolve(ode,y(x)); else sol := :-dsolve([ode,ic],y(x)); fi; end proc; export get_sol::static:=proc(_self,$) return sol; end proc; end module;
And now use it as follows
o:=Object(ODE, diff(y(x),x)+y(x)=sin(x), y(x), y(0)=0); o:-dsolve(): o:-get_sol(); #y(x) = -1/2*cos(x) + 1/2*sin(x) + 1/2*exp(-x) o:=Object(ODE, diff(y(x),x)+y(x)=sin(x), y(x)); o:-dsolve(): o:-get_sol(); #y(x) = -1/2*cos(x) + 1/2*sin(x) + exp(-x)*_C1
In the child class you want to extend from the parent class, add option object(ParentName);
Here is an example
restart; module ODE() option object; local ode; export set_ode::static:=proc(_self,ode,$) _self:-ode :=ode; end proc; export get_ode::static:=proc(_self,$) return _self:-ode; end proc; end module; #create class/module which extends the above module second_order_ode() option object(ODE); export get_ode_order::static:=proc(_self,$) return 2; end proc; end module;
In the above second_order_ode inherts all local variables and functions in the ODE class and
adds new proc. Use as follows
o:=Object(second_order_ode); #create an object instance o:-set_ode(diff(y(x),x)=sin(x)); o:-get_ode(); o:-get_ode_order();
Note that the child class can not have its own variable with the same name as the parent
class. This is limitation. in C++ for example, local variables in extended class overrides the
same named variable in the parent class.
Even if the variable have different type, Maple will not allow overriding. For example, this will fail
restart; module ODE() option object; local ode; local id::integer; export set_ode::static:=proc(_self,ode,$) print("Enter ode::set_ode"); _self:-ode :=ode; end proc; export get_ode::static:=proc(_self,$) return _self:-ode; end proc; end module; module second_order_ode() option object(ODE); local id::string; export get_ode_order::static:=proc(_self,$) return 2; end proc; end module; Error, (in second_order_ode) local `id` is declared more than once
There might be a way to handle this, i.e. to somehow exlicitly tell Maple to override parant class proc or variable name in the child. I do not know now. The above is using Maple 2021.1
This is called polymorphism in OOP. This is a base class animal_class which has
make_sound method. This method acts just as a place holder (interface), which the
extending class must extends (override) with an actual implementation.
The class is extended to make cat class and implementation is made.
module animal_class() option object; export make_sound::static:=proc(_self,$) error("Not implemented, must be overriden"); end proc; end module; %--------------- module cat_class() option object(animal_class); make_sound::static:=proc(_self,$) #note. DO NOT USE export print("mewooo"); end proc; end module;
And now
o:=Object(animal_class); o:-make_sound(); Error, (in anonymous procedure) Not implemented, must be overriden
The above is by design. As the animal_class is meant to be extended to be usable.
my_cat:=Object(cat_class); my_cat:-make_sound(); "mewooo"
So the base class can have number of methods, which are all meant to be be have its implementation provided by an extending class. Each class which extends the base class must provide implementation.
Once a base class is extended, all methods in the base class become part of the extending class. So to call a base class method just use same way as if calling any other method in the extending class itself.
Here is an example.
module person() option object; local age::integer:=100; export ModuleCopy::static:= proc(_self,proto::person, age::integer,$) _self:-age := age; end proc; local base_class_method::static:=proc(_self,$) print("In base class method..."); end proc; end module; #---- extend the above class module young_person() option object(person); export process::static:=proc(_self,$) print("In young_person process"); _self:-base_class_method(); end proc; end module;
Here is an example of usage
o:=Object(young_person,20); o:-process() "In young_person process" "In base class method..."
The above in Maple 2023.1
A Maple Object can be used a record type in other languags, such as Ada or Pascal. This example shows how to define a local type inside a proc and use it as record.
restart; foo:=proc(the_name::string,the_age::integer)::person_type; local module person_type() #this acts as a record type option object; export the_name::string; export the_age::integer; end module; local person::person_type:=Object(person_type); person:-the_name:=the_name; person:-the_age:=the_age; return person; end proc; o:=foo("joe doe",100); o:-the_name; "joe doe" o:-the_age; 100
In the above person is local variable of type person_type. In the above example, the local
variable was returned back to user. But this is just an example. One can declare such
variables and just use them internally inside the proc only. This method helps one organize
related variables into one record/struct like type. The type can also be made global if
needed.
Suppose we have list L1 of objects and we want to copy this list to another list, say
L2. If we just do L2 = L1 then this will not make an actual copy as any changes
to L1 are still reflected in L2. The above only copies the reference to the same
list.
To make a physical copy, need to use the copy command as follows
restart; module person_type() option object; export the_name::string:="doe joe"; export the_age::integer:=100; end module; L1:=[]; for N from 1 to 5 do o:=Object(person_type); o:-the_name:=convert(N,string); L1:= [ op(L1), o]; od: L2:=map(Z->copy(Z),L1):
Now making any changes to L1 will not affect L2. If we just did L2 = L1 then both will
share same content which is not what we wanted.
This is a basic example of using OOP in Maple to implement an ode solver. There is a base module called ode_base_class (I will be using class instead of module, as this is more common in OOP)
This program will for now support first order and second order ode’s only.
The base ode class will contain the basic operations and data which is common to both first order and second order ode’s.
Next, we will make a first order ode class, and second order ode class. Both of these will extend the base ode class.
Next, we will have more ode classes. For example, for first order ode, there will be linear first order ode class, and separable first order ode class, and Bernoulli ode class and so on. Each one these classess will extend the first order ode class which in turn extends the base ode class.
Same for second order ode’s. There will be second order constant coefficients ode class, and second order variable coefficient ode class and so on. Each one of these classes will extend the second order ode class which in turn extends the base ode class.
Let base_ode_class be the base of an ode class which will contain all the neccessary
methods that are generic and applicable to an ode of any order and type.
These can be the ode expression itself, its order and any other data and operations which applicable to any ode type.
Now we want to create a first order class. This will have its own operations and private variables that are specific and make sense only to any first order ode.
Then these will be the first order separable ode class, which has methods that implement solving the ode using separable method and has other methods which makes sense only for first order separable ode. The following diagram is partial illustration of the is-A relation among possible classes.
First we define the base ode class and define private variables and method that are common to any ode type.
restart; module base_ode_class() option object; local the_ode; local the_order::integer; export get_ode::static:=proc(_self,$) RETURN(_self:-the_ode); end proc; export get_order::static:=proc(_self,$) RETURN(_self:-the_order); end proc; end module;
Note that the base ode class does not have constructor. Since it is meant to be extended only.
The following is the first order ode class.
module first_order_ode_class() option object(base_ode_class); local initial_conditions; local solution_found; export get_IC::static:=proc(_self,$) RETURN(_self:-initial_conditions); end proc; export get_solution::static:=proc(_self,$) RETURN(_self:-solution_found); end proc; export verify_solution::static:=proc(_self,$)::truefalse; #code to verify if the solution found is valid or not #using odetest() end proc; end module;
The following is the first order separable ode class which extends the above first order ode class.
module first_order_separable_ode_class() option object(first_order_ode_class); local f,g; #ode of form y'=f(x)*g(y) export ModuleCopy::static:= proc(_self,proto::first_order_separable_ode_class, ode, IC,$) _self:-the_ode := ode; _self:-initial_conditions := IC; end proc; export dsolve::static:=proc(_self,$) #solve the ode for now we will use Maple but in my code #I have my own solver ofcourse. _self:-solution_found:= :-dsolve([_self:-the_ode, _self:-initial_conditions]); end proc; end module;
In the above, when we create a instance of first_order_separable_ode_class then it now
have the whole chain of classes into one. i.e. first order separable class extending the first
order class which in turn extends the base ode class. For example
o:=Object(first_order_separable_ode_class,diff(y(x),x)=3*sin(x)*y(x),y(0)=1) o:-dsolve(); o:-get_solution() #y(x) = exp(3)*exp(-3*cos(x)) o:-get_ode() #diff(y(x), x) = 3*sin(x)*y(x)
The above calls will all work, even though the first order separable class has no
method in it called get_ode but it extends a class which does, hence it works as
it.
Now we will do the same for second order ode’s.
Advantage of this design, is that methods in base classes that are being extended can be reused by any class which extends them. Only methods that applies and specific to the lower level classes need to be implemented.
As we add more specific solvers, we just have to extend the base classes and the new solvers
just need to implement its own specific dsolve and any specific methods and data that it
needs itself.
Ofcourse in practice the above design is not complete as is. The user should not have to specify which class to instantiate, as user does not care what the ode type or class it is. They just want to to do
o:=Object(ode_class,diff(y(x),x)=3*sin(x)*y(x),y(0)=1) o:-dsolve(); o:-get_solution()
To solve this problem we have to make a factory method which is called to make the correct instance of the class and return that to the user. The factory method figures out the type of ode and it creates the correct instance of the correct class and returns that. So the call will become
o := make_ode_object( diff(y(x),x)=3*sin(x)*y(x), y(0)=1) o:-dsolve(); o:-get_solution()
The function make_ode_object above is the main interface the user will call to make an ode
object.
This will be explained next with examples. One possibility is to make the factory function a global function or better, a public method in a utility module. For now, it is given here as stadalone function for illustration. The user calls this method to make an object of the correct instance of the ode. Here is complete implementation of all the above including the factory method.
#factory method. Makes objects for users make_ode_object:=proc(ode::`=`,func::function(name)) local x,y,the_order; y:=op(0,func); x:=op(1,func); the_order := PDEtools:-difforder(ode,x); if the_order=1 then RETURN(first_order_ode_class:-make_ode_object(ode,func)); elif the_order=2 then #RETURN(second_order_ode_class:-make_ode_object(ode,func)); #implement later NULL; else error "Only first and second order ode's are currently supported"; fi; end proc: ######################### module base_ode_class() option object; local the_ode; local the_order::integer; #methods bound to the object export get_ode::static:=proc(_self,$) RETURN(_self:-the_ode); end proc; export get_order::static:=proc(_self,$) RETURN(_self:-the_order); end proc; end module: ################## module first_order_ode_class() option object(base_ode_class); local initial_conditions; local solution_found; #public factory method not bound to the object. export make_ode_object:=proc(ode::`=`,func::function(name)) local x,y,ode_type::string; y:=op(0,func); x:=op(1,func); ode_type:="separable"; #code here which determined first order ode type if ode_type="separable" then RETURN( Object(first_order_separable_ode_class,ode,func)); elif ode_type="linear" then RETURN( Object(first_order_linear_ode_class,ode,func)); fi; #more ode types added here end proc; #methods bound to the object export get_IC::static:=proc(_self,$) RETURN(_self:-initial_conditions); end proc; export get_solution::static:=proc(_self,$) RETURN(_self:-solution_found); end proc; export verify_solution::static:=proc(_self,$)::truefalse; #code to verify if the solution found is valid or not #using odetest() end proc; end module: ################## module first_order_separable_ode_class() option object(first_order_ode_class); local f,g; #ode of form y'=f(x)*g(y) export ModuleCopy::static:= proc(_self,proto::first_order_separable_ode_class,ode,func::function(name),$) _self:-the_ode := ode; end proc; export dsolve::static:=proc(_self,$) #print("Enter first_order_separable_ode_class:-dsolve"); #solve the ode for now we will use Maple but in my code #I have my own solver ofcourse. _self:-solution_found:= :-dsolve(_self:-the_ode); NULL; end proc; end module:
It is used as follows
o:=make_ode_object(diff(y(x),x)=sin(x)*y(x),y(x)); o:-dsolve(); o:-get_solution(); y(x) = c__1 exp(-cos(x))
Here, I will start making a complete small OOP ode solver in Maple.
At each step more classes are added and enhanced until we get a fully working small ode solver based on OOP design that solves a first and second order ode, this is to show how it all works. Another solvers can be added later by simply extending the base class.
The base class is called Base_ode_class. There will be Second_order_ode_class and
First_order_ode_class and these classes extend the base ode class. We can later add
higher_order_ode_class.
Next, there are different classes which extend these. There is First_order_linear_ode_class
and First_order_separable_ode_class and so on, and these extend the
First_order_ode_class.
For example, if a user wanted to solve a first order ode which happend to be say
separable, then object of class First_order_separable_ode_class will be created and
used.
Since the user does not know and should not know what object to create, then the factory class will be used. The factory class is what the user initially calls to make the ode object.
It is the factory class which determines which type of class to instantiate based on the ode given.
The factory class is singleton (standard module in Maple, not of type object),
which has the make_ode method which is called by the user. This method parses
the ode and determines its order and then based on the order determine which
subclass to use, and then instantiate this and returns the correct object to the user to
use.
This object will have the dsolve method and other methods the user can use on the
object.
The make_ode method in the factory module accepts only the ode itself the function such as
\(y(x)\). A typical use is given below
ode := ODE_factoy_class:-make_ode( diff(y(x),x)=sin(x), y(x) ); ode:-set_IC(....); ode:-set_hint("the hint"); ..... ode:-dsolve(); #solves the ode ode:-is_solved(); #checks if ode was successfully solved ode:-verify_sol(); #verifies the solution using maple odetest() ode:-is_sol_verified(); #asks if solution is verified ode:-get_number_of_sol(); #returns number of solutions, 1 or 2 etc... ode:-get_sol(); #returns the solutions found in list #and more method ...
Examples at the end will show how all the above works on actual odes’s.
The initial call to make an ode does not have initial conditions, or hint and any other parameters. This is so to keep the call simple. As the factory method only makes the concrete ode object.
Additional methods are then used to add more information if needed by using the returned object itself, such as initial conditions, and hint and so on before calling the dsolve method on the object.
Here is a very basic setup which include the base ode class and extended to first order two subclasses for now.
restart; ODE_factory_class :=module() #notice, normal module. No option object. export make_ode:=proc(ode::`=`,func::function(name),$) local dep_variables_found::list,item; local y::symbol; local x::symbol; local ode_order::integer; if nops(func)<>1 then error("Parsing error, dependent variable must contain one argument, found ", func); fi; y:=op(0,func); x:=op(1,func); if not has(ode,y) then error ("Supplied ode ",ode," has no ",y); fi; if not has(ode,x) then error ("Supplied ode ",ode," has no ",x); fi; if not has(ode,func) then error ("Supplied ode ",ode," has no ",func); fi; ode_order := PDEtools:-difforder(ode,x); #this will check that the dependent variable will show with #SAME argument in the ode. i.e. if y(x) and y(t) show up in same ode, it #will throw exception, which is what we want. try dep_variables_found := PDEtools:-Library:-GetDepVars([y],ode); catch: error lastexception; end try; #now go over dep_variables_found and check the independent #variable is same as x i.e. ode can be y'(z)+y(z)=0 but function is y(x). for item in dep_variables_found do if not type(item,function) then error("Parsing error. Expected ",func," found ",item," in ode"); else if op(1,item) <> x then error("Parsing error. Expected ",func," found ",item," in ode"); fi; fi; od; #now go over all indents in ode and check that y only shows as y(x) and not as just y #as the PDEtools:-Library:-GetDepVars([_self:-y],ode) code above does not detect this. #i.e. it does not check y'(x)+y=0 if numelems(indets(ode,identical(y))) > 0 then error("Parsing error, Can not have ",y," with no argument inside ",ode); fi; if ode_order=1 then RETURN(make_first_order_ode(ode,y,x)); elif ode_order=2 then RETURN(make_second_order_ode(ode,y,x)); else RETURN(make_higher_order_ode(ode,y,x)); fi; end proc; local make_first_order_ode:=proc(ode::`=`,y::symbol,x::symbol) #decide on what specific type the ode is, and make instant of it if first_order_ode_quadrature_class:-is_quadrature(ode,y,x) then RETURN(Object(first_order_ode_quadrature_class,ode,y,x)); elif first_order_ode_linear_class:-is_linear(ode,y,x) then RETURN(Object(first_order_ode_linear_class,ode,y,x)); fi; #and so on end proc; local make_second_order_ode:=proc(ode::`=`,y::symbol,x::symbol) #decide on what specific type the ode is, and make instant of it #same as for first order end proc; end module; #------------------------- module solution_class() option object; local the_solution; local is_verified_solution::truefalse:=false; local is_implicit_solution::truefalse:=false; export ModuleCopy::static:= proc(_self,proto::solution_class,the_solution::`=`,is_implicit_solution::truefalse,$) _self:-the_solution:=the_solution; _self:-is_implicit_solution:=is_implicit_solution; end proc; export get_solution::static:=proc(_self,$) RETURN(_self:-the_solution); end proc; export is_verified::static:=proc(_self,$) RETURN(_self:-is_verified_solution); end proc; export is_implicit::static:=proc(_self,$) RETURN(_self:-is_implicit_solution); end proc; export is_explicit::static:=proc(_self,$) RETURN(not(_self:-is_implicit_solution)); end proc; export verify_solution::static:= overload( [ proc(_self, ode::`=`,$) option overload; local stat; stat:= odetest(_self:-the_solution,ode); if stat=0 then _self:-is_verified_solution:=true; else if simplify(stat)=0 then _self:-is_verified_solution:=true; else _self:-is_verified_solution:=false; fi; fi; end, proc(_self, ode::`=`,IC::list,$) option overload; local stat; stat:= odetest([_self:-the_solution,IC],ode); if stat=[0,0] then _self:-is_verified_solution:=true; else if simplify(stat)=[0,0] then _self:-is_verified_solution:=true; else _self:-is_verified_solution:=false; fi; fi; end ]); end module: #------------------------- module ODE_base_class() option object; local y::symbol; local x::symbol; local func::function(name); #y(x) local ode::`=`; local ode_order::posint; local IC::list:=[]; local parsed_IC::list:=[]; local the_hint::string:=""; local solutions_found::list(solution_class):=[]; #exported getters methods export get_ode::static:=proc(_self,$) RETURN(_self:-ode); end proc; export get_x::static:=proc(_self,$) RETURN(_self:-x); end proc; export get_y::static:=proc(_self,$) RETURN(_self:-y); end proc; export get_ode_order::static:=proc(_self,$) RETURN(_self:-ode_order); end proc; export get_IC::static:=proc(_self,$) RETURN(_self:-IC); end proc; export get_parsed_IC::static:=proc(_self,$) RETURN(_self:-parsed_IC); end proc; export get_sol::static:=proc(_self,$) local L:=Array(1..0): local sol; for sol in _self:-solutions_found do L ,= sol:-get_solution(); od; RETURN(convert(L,list)); end proc; #exported setters methods export set_hint::static:=proc(_self,hint::string,$) #add code to check if hint is valid _self:-the_hint:=hint; end proc; end module; #------------------------- module first_order_ode_quadrature_class() option object(ODE_base_class); local f,g; #ode of form y'=f(x)*g(y) #this method is not an object method. It is part of the module but does #not have _self. It is called by the factory class to find if the ode #is of this type first export is_quadrature:=proc(ode::`=`,y::symbol,x::symbol)::truefalse; RETURN(true); #for now end proc; export ModuleCopy::static:= proc(_self,proto::first_order_ode_quadrature_class,ode::`=`,y::symbol,x::symbol,$) _self:-ode := ode; _self:-y := y; _self:-x := x; _self:-func := _self:-y(_self:-x); _self:-ode_order :=1; end proc; export dsolve::static:=proc(_self,$) local sol,o; #print("Enter first_order_ode_quadrature_class:-dsolve"); #solve the ode for now we will use Maple but in my code #I have my own solver ofcourse. sol:= :-dsolve(_self:-ode,_self:-func); o:=Object(solution_class,sol,false); _self:-solutions_found:= [o]; NULL; end proc; end module: #------------------------- module first_order_ode_linear_class() option object(ODE_base_class); local f,g; #ode of form y'=f(x)*g(y) #this method is not an object method. It is part of the module but does #not have _self. It is called by the factory class to find if the ode #is of this type first export is_linear:=proc(ode::`=`,y::symbol,x::symbol)::truefalse; RETURN(true); #for now end proc; export ModuleCopy::static:= proc(_self,proto::first_order_ode_linear_class,ode::`=`,y::symbol,x::symbol,$) _self:-ode := ode; _self:-y := y; _self:-x := x; _self:-func := _self:-y(_self:-x); _self:-ode_order :=1; end proc; export dsolve::static:=proc(_self,$) local sol,o; sol:= :-dsolve(_self:-ode,_self:-func); o:=Object(solution_class,sol,false); _self:-solutions_found[1]:= [o]: end proc: end module:
Example usage is
o:=ODE_factory_class:-make_ode(diff(y(x),x)=x,y(x)) o:-get_ode() d --- y(x) = x dx o:-dsolve(); o:-get_sol() [ 1 2 ] [y(x) = - x + c__1] [ 2 ]
The commnand DEtools:-odeadvisor(ode,y(x)); gives the ode type names.
See https://maplesoft.com/support/help/maple/view.aspx?path=DEtools/odeadvisor for known names of ode type in Maple.
But we can also tell it to check if the ode is of specific type using
DEtools:-odeadvisor(ode,y(x),[name]);. For an example
ode:=x^2+3*x*diff(y(x),x)=y(x)^4+2*y(x); DEtools:-odeadvisor(ode,y(x),[Chini]) [_Chini]
If the name we have given is not the type of the ode, Maple will return [NONE].
For an example
ode:=x^2+3*x*diff(y(x),x)=y(x)^4+2*y(x); DEtools:-odeadvisor(ode,y(x),[separable]) [NONE]
To find what methods dsolve uses do
`dsolve/methods`[1] [quadrature, linear, Bernoulli, separable, inverse_linear, homogeneous, Chini, lin_sym, exact, Abel, pot_sym]
The above gives methods for first order ODEs. More are given by
`dsolve/methods`[1,'semiclass'] [Riccati, inverse_Riccati, equivalent_to_Abel, linearizable, linearizable_by_differentiation] `dsolve/methods`[1,'high_degree'] [WeierstrassP, WeierstrassPPrime, JacobiSN, linearizable_by_differentiation, missing, dAlembert, homogeneous_B, sym_implicit] `dsolve/methods`[1,"development"] [linearizable_by_differentiation, linearizable, con_sym, WeierstrassP, WeierstrassPPrime, equivalent_to_Abel, Abel_AIR, special, Riccati_symmetries] `dsolve/methods`[1,extra] [inverse_Riccati, Abel_AIL, `sym_pat/[F(x)*G(y),0]`, `sym_pat/[F(x),G(x)]`, `sym_pat/[F(x),G(y)]`, `sym_pat/[F(x)+G(y),0]`, `sym_pat/[F(x),G(x)*y+H(x)]`, sym_pat, exp_sym]
For example given a first order ode, we can ask dsolve to solve it using one of these methods as follows
restart; ode:=diff(y(x),x)=sqrt( sqrt( y(x) )* sin(x))/sqrt(y(x)) dsolve(ode,y(x),[`exact`]) 4/5*y(x)^(5/4) + Intat(-sqrt(sqrt(y(x))*sin(_a))/y(x)^(1/4), _a = x) + c__1 = 0
We can ask it to use symmetry with specific pattern as follows
restart; ode:=diff(y(x),x)=sqrt( sqrt( y(x) )* sin(x))/sqrt(y(x)) dsolve(ode,y(x),[`sym_pat/[F(x),G(x)*y+H(x)]`]) c__1 + 4/5*y(x)^(3/2)*sqrt(sin(x))/sqrt(sqrt(y(x))*sin(x)) - Int(sqrt(sin(x)), x) = 0 #or using short cut, same thing can be done as follows dsolve(ode,y(x),Lie) c__1 + 4/5*y(x)^(3/2)*sqrt(sin(x))/sqrt(sqrt(y(x))*sin(x)) - Int(sqrt(sin(x)), x) = 0 #another short cut for Lie is dsolve(ode,[`sym_pat`]) c__1 + 4/5*y(x)^(3/2)*sqrt(sin(x))/sqrt(sqrt(y(x))*sin(x)) - Int(sqrt(sin(x)), x) = 0
To find all methods, and for higher order ode, we first use indices as follows
indices(`dsolve/methods`) [high, linear_nonhomogeneous], [1], [1, high_degree], [3, linear_homogeneous], [2, "linear_homogeneous all methods"], [2, "linear_homogeneous other"], [2, linear_homogeneous], [2, "special_functions"], [3, linear_nonhomogeneous], [2, "linear_homogeneous in Normal Form"], [2, "development"], [2, "hypergeometric"], [1, extra], [high, linear_homogeneous], [high, nonlinear], [1, "special"], [1, "development"], [2, nonlinear], [high, "development"], [1, semiclass], [3, nonlinear], [2, linear_nonhomogeneous], [2, "linear_homogeneous as given"], [3, "development"]
Using the above, to find all methods for second order linear_homogeneous ode, the command is
`dsolve/methods`[2, "linear_homogeneous all methods"] [quadrature, const_coeffs, Euler, linear_1, `linear/missing_y`, Kovacic, special_functions, to_const_coeffs, exact_linear, sym_1, Mathieu, MeijerG, Heun, HeunG, HeunC, HeunB, HeunD, HeunT, mu_xy, equivalent_to_Bessel, to_Riccati, Bessel, elliptic, Legendre, Whittaker, Kummer, cylindrical, hypergeometric, hypergeom1, hypergeom2, Riemann, RNF, hypergeometricsols, rationalize_lode, with_periodic_functions] `dsolve/methods`[2, "linear_homogeneous other"] [exact_linear, sym_1, to_const_coeffs, mu_xy, equivalent_to_Bessel, to_Riccati, with_periodic_functions]
For example, given the ode \(y''+3 y'+ y=0\), we can now do
dsolve(ode,y(x)) y(x) = c__1*exp(1/2*(sqrt(5) - 3)*x) + c__2*exp(-1/2*(3 + sqrt(5))*x) dsolve(ode,y(x),[`const_coeffs`]) y(x) = c__1*exp(1/2*(sqrt(5) - 3)*x) + c__2*exp(-1/2*(3 + sqrt(5))*x) dsolve(ode,y(x),[`Kovacic`]) y(x) = c__1*exp(1/2*(sqrt(5) - 3)*x) + c__2*exp(-1/2*(3 + sqrt(5))*x)
Not all methods ofcourse will work, as it depends on the ode type.
This function below lists all methods
ind:=indices(`dsolve/methods`); for item in ind do cat("`dsolve/methods`",String(item)); eval(parse(%)) od;
Which gives
"`dsolve/methods`[high, linear_nonhomogeneous]" [quadrature, fully_exact_linear, linear_nonhomogeneous_[0,1], exact_linear_nonhomogeneous, linear, exp_reduce] "`dsolve/methods`[1]" [quadrature, linear, Bernoulli, separable, inverse_linear, homogeneous, Chini, lin_sym, exact, Abel, pot_sym] "`dsolve/methods`[1, high_degree]" [WeierstrassP, WeierstrassPPrime, JacobiSN, linearizable_by_differentiation, missing, dAlembert, homogeneous_B, sym_implicit] "`dsolve/methods`[3, linear_homogeneous]" [quadrature, const_coeffs, Euler, fully_exact_linear, to_const_coeffs, linear, exp_reduce, exact_linear, with_periodic_functions] "`dsolve/methods`[2, "linear_homogeneous all methods"]" [quadrature, const_coeffs, Euler, linear_1, linear/missing_y, Kovacic, special_functions, to_const_coeffs, exact_linear, sym_1, Mathieu, MeijerG, Heun, HeunG, HeunC, HeunB, HeunD, HeunT, mu_xy, equivalent_to_Bessel, to_Riccati, Bessel, elliptic, Legendre, Whittaker, Kummer, cylindrical, hypergeometric, hypergeom1, hypergeom2, Riemann, RNF, hypergeometricsols, rationalize_lode, with_periodic_functions] "`dsolve/methods`[2, "linear_homogeneous other"]" [exact_linear, sym_1, to_const_coeffs, mu_xy, equivalent_to_Bessel, to_Riccati, with_periodic_functions] "`dsolve/methods`[2, linear_homogeneous]" [linear_homogeneous] "`dsolve/methods`[2, "special_functions"]" [Bessel, elliptic, Legendre, Kummer, Whittaker, hypergeometric, Mathieu] "`dsolve/methods`[3, linear_nonhomogeneous]" [quadrature, fully_exact_linear, linear_nonhomogeneous_[0,1], exact_linear_nonhomogeneous, linear, exp_reduce] "`dsolve/methods`[2, "linear_homogeneous in Normal Form"]" [linear_1] "`dsolve/methods`[2, "development"]" [mu_xyp, mu_xyp2, mu_formal, mu_heuristic, exp_reduce, linear, Bessel2, Whittaker_old, nonlinear_homogeneous, exact_linear_nonhomogeneous, mu_y1, mu_x_y1, mu_y_y1, mu_poly_yn, exp_sym, sym_pat, sym_8] "`dsolve/methods`[2, "hypergeometric"]" [hypergeom1, hyper3] "`dsolve/methods`[1, extra]" [inverse_Riccati, Abel_AIL, sym_pat/[F(x)*G(y),0], sym_pat/[F(x),G(x)], sym_pat/[F(x),G(y)], sym_pat/[F(x)+G(y),0], sym_pat/[F(x),G(x)*y+H(x)], sym_pat, exp_sym] "`dsolve/methods`[high, linear_homogeneous]" [quadrature, const_coeffs, Euler, fully_exact_linear, to_const_coeffs, linear, exp_reduce, exact_linear, with_periodic_functions] "`dsolve/methods`[high, nonlinear]" [linearizable_by_differentiation, linearizable, reducible, exact_nonlinear, missing, mu_formal, lin_sym] "`dsolve/methods`[1, "special"]" [80, 81] "`dsolve/methods`[1, "development"]" [linearizable_by_differentiation, linearizable, con_sym, WeierstrassP, WeierstrassPPrime, equivalent_to_Abel, Abel_AIR, special, Riccati_symmetries] "`dsolve/methods`[2, nonlinear]" [Liouville, WeierstrassP, JacobiSN, linearizable, linearizable_by_differentiation, mu_xy_2, missing, mu_xyp2_dynamical_symmetries_fully_reducible, mu_xyp_singularcases, sym_1, exact_nonlinear, reducible, lin_sym, S-function, mu_xyp_generalcase, mu_xyp2_dynamical_symmetries_not_fully_reducible] "`dsolve/methods`[high, "development"]" [k25, RNF, mu_heuristic, MeijerG, nonlinear_homogeneous, mu_poly_yn, exp_sym] "`dsolve/methods`[1, semiclass]" [Riccati, inverse_Riccati, equivalent_to_Abel, linearizable, linearizable_by_differentiation] "`dsolve/methods`[3, nonlinear]" [linearizable_by_differentiation, linearizable, missing, exact_nonlinear, reducible, mu_formal, lin_sym] "`dsolve/methods`[2, linear_nonhomogeneous]" [quadrature, fully_exact_linear, linear_nonhomogeneous_[0,1], linear_nonhomogeneous_[0,F(x)], linear_nonhomogeneous] "`dsolve/methods`[2, "linear_homogeneous as given"]" [quadrature, const_coeffs, Euler, linear_1, linear/missing_y, Kovacic, RNF, special_functions, MeijerG, Heun, hypergeometricsols, rationalize_lode] "`dsolve/methods`[3, "development"]" [k25, RNF, mu_heuristic, linear_patterns, MeijerG, nonlinear_homogeneous, mu_y2, mu_poly_yn, exp_sym, pFq, 3F2, 2F2, 1F2, 0F2]
restart; ode:=diff(y(x),x)+y(x)^2*sin(x)-2*sin(x)/cos(x)^2 = 0; yp:=DETools:-particularsol(ode);
To step into the code, do
restart; ode:=diff(y(x),x)+y(x)^2*sin(x)-2*sin(x)/cos(x)^2 = 0; stopat(`DEtools/particularsol`); DETools:-particularsol(ode);
To print it do
print(`DEtools/particularsol`);
Use the output=basis option
ode:=diff(y(x),x$2)-x*diff(y(x),x)-x*y(x)=0; dsolve(ode,output=basis);
To solve \[ y''-3y'+2y=10 e^{5 x} \] with \(y(0)=1,y'(0)=5\) do
eq1:= diff(y(x),x$2)-3*diff(y(x),x)+2*y(x)=10*exp(5*x); dsolve({eq1,y(0)=1,D(y)(0)=5},y(x)); Methods for second order ODEs: Trying to isolate the derivative d^2y/dx^2... Successful isolation of d^2y/dx^2 --- Trying classification methods --- trying a quadrature trying high order exact linear fully integrable trying differential order: 2; linear nonhomogeneous with symmetry [0,1] trying a double symmetry of the form [xi=0, eta=F(x)] <- double symmetry of the form [xi=0, eta=F(x)] successful ....
The above can also be written using D@@ notation, like this
eq:= (D@@2)(y)(x) - 3*D(y)(x) +2*y(x) = 10*exp(5*x); IC := y(0)=1,D(y)(0)=5; dsolve({eq,IC},y(x));
use odetest and check if it gives zero.
eq1:= diff(diff(y(x),x),x)-3*diff(y(x),x)+2*y(x)=10*exp(5*x); ans:=dsolve({eq1,IC},y(x)); odetest(ans,eq1); 0
Maple can classify the ODE.
eq1:= diff(y(x),x$2)-3*diff(y(x),x)+2*y(x)=10*exp(5*x); R0 := DEtools['odeadvisor'](eq1,y(x)); R0 := [[_2nd_order, _with_linear_symmetries]]
To get help on this type of ODE, do
DEtools['odeadvisor'](eq1,'help');
Use with(DEtools);
restart; eq1:= diff(y(x),x$2)-3*diff(y(x),x)+2*y(x)=10*exp(5*x); DEtools[DEplot](eq1,y(x),x=-2..5, [ [y(0)=0, D(y)(0)=0]], y=-3..3,linecolor=red);
To get a better plot, change the stepsize and independent variable range
restart; eq1:= diff(y(x),x$2)-3*diff(y(x),x)+2*y(x)=10*exp(5*x); DEtools[DEplot](eq1,y(x),x=-1..1,[[y(0)=0,D(y)(0)=0]],y=-3..3,stepsize=0.001,linecolor=red);
From: Robert Israel (israel@math.ubc.ca) Subject: Re: given precision in Maple Newsgroups: comp.soft-sys.math.maple Date: 2003-07-16 20:19:06 PST Set Digits:= n and all calculations from this point will be done with n digits. Mathematical functions will be correct to n digits as well (to the extent this is practical). If you want high-accuracy numerical ODE solutions, on the other hand, it's not so simple. I think the best way is using the taylorseries method. For example, consider the problem y' = y^2, y(1) = 1, where the exact solution y = 1/(2-x) has y(1.9) = 10. > Digits:= 30: sol:= dsolve({D(y)(x)=y(x)^2, y(1) = 1}, y(x), numeric, method=taylorseries, abserr=1e-25): sol(1.9); [x = 1.9, y(x) = 9.99999999999999999999999797691] > 10 - eval(y(x),%); -23 0.202309 10 The other methods (in particular the default rkf45) do not give results anywhere near this good.
How to obtain the ODE in canonical coordinates \(S(R)\) for an ode using Lie symmetry?
Here is a function I wrote with the help of Maple docs, that does that. For now, I am using first order ode. It takes the ode and the variables (coordinates) to use (\(R,S\)) is what normally used, and returns back the ODE in canonical coordinates. The ode should always be a quadrature for first order.
get_ODE_in_canonical:=proc(ode::`=`,y::symbol,x::symbol,S::symbol,R::symbol) local infinitesimals,tr,itr,ODE; infinitesimals := [DEtools:-symgen(ode)]; if nops(infinitesimals)=0 then RETURN(FAIL); fi; infinitesimals:=infinitesimals[1]; tr := DEtools:-canoni(infinitesimals,y(x),S(R)); itr := op(1,[solve(tr,{x,y(x)})]); ODE:= PDEtools:-dchange(itr,ode,[R,S(R)],simplify); ODE:= op(solve(ODE,{diff(S(R),R)})); RETURN(ODE); end proc:
Example
ODE := diff(y(x),x) = 1/x^2+y(x); get_ODE_in_canonical(ODE,y,x,S,R) #diff(S(R), R) = 1/(exp(R)*R^2) ODE := diff(y(x),x) = x^2-2*x*y(x)+y(x)^2; get_ODE_in_canonical(ODE,y,x,S,R) #diff(S(R), R) = 1
The following function is small change of the above. It accepts your choice of \(\xi ,\eta \) to use instead of the first one obtained by Maple. This makes it possible to experiment with different infinitesimals. Note that different infinitesimals will/can result in different canonical ODE, but the final solution for \(y(x)\) of course will always be the same after applying the transformation back to \(x,y\) coordinates.
restart; get_ODE_in_canonical_v2:=proc(ode::`=`,y::symbol,x::symbol,S::symbol,R::symbol,xi,eta) local infinitesimals:=[xi,eta],tr,itr,ODE; tr := DEtools:-canoni(infinitesimals,y(x),S(R)); itr := op(1,[solve(tr,{x,y(x)})]); ODE:= PDEtools:-dchange(itr,ode,[R,S(R)],simplify); ODE:= op(solve(ODE,{diff(S(R),R)})); RETURN(simplify(ODE,symbolic)); end proc:
And now
ODE:=diff(y(x),x) = x^2-2*x*y(x)+y(x)^2; get_ODE_in_canonical_v2(ODE,y,x,S,R,1,1) #diff(S(R), R) = 1/(R^2 - 1)
Here is an example.
ode:=diff(y(x),x)=-y(x)^2/(exp(x)-y(x)) DEtools:-symgen(ode,y(x),HINT=[c1*x+c2*y+c3,c4*x+c5*y+c6]); [_xi = 1, _eta = y]
Or
ode:=diff(y(x),x)=-y(x)^2/(exp(x)-y(x)) DEtools:-symgen(ode,y(x),HINT=[g(x),f(x)*y]); [_xi = 1, _eta = y]
So not use y(x) but only y in the hint.
The following function takes in a single ode and parses it to verify it is valid. It returns back 3 values. The dependent variable, the indepenent variable and the ode order.
#added 12/8/2022. interface(warnlevel=4); kernelopts('assertlevel'=2): parse_single_ode:=proc(ode::`=`)::symbol,symbol,integer; #parses single ode. returns back dep_var,indep_var,ode_order #throws exception when it detects parsing errors local func; local y,x; local dep_variables_found::list,item; local the_order; func:=PDEtools:-Library:-GetDepVars("not given",ode,onlydifferentiatedfunctions=true); if nops(func)=0 then error ("not differential equation ",ode); fi; func := func[1]; if nops(func)<>1 then error("Parsing error, dependent variable must contain one argument, found ", func); fi; y:=op(0,func); x:=op(1,func); #basic verification if not has(ode,y) then error ("Supplied ode ",ode," has no ",y); fi; if not has(ode,x) then error ("Supplied ode ",ode," has no ",x); fi; if not has(ode,func) then error ("Supplied ode ",ode," has no ",func); fi; the_order := PDEtools:-difforder(ode,x); #if the_order=0 then # error ("No derivative found in ",ode,". Input is not differential equation"); #fi; #note that the following call will also return y(x) if the input is not an ode #this will check that the dependent variable will show with SAME argument in the ode #i.e. if y(x) and y(t) show up in same ode, it will throw exception, which is what #we want. try dep_variables_found := PDEtools:-Library:-GetDepVars([y],ode); catch: error lastexception; end try; #now go over dep_variables_found and check the independent variable is same as x #i.e. ode can be y'(z)+y(z)=0 but function is y(x). for item in dep_variables_found do if not type(item,function) then error("Parsing error. Expected ",func," found ",item," in ode"); else if op(1,item) <> x then error("Parsing error. Expected ",func," found ",item," in ode"); fi; fi; od; #now go over all indents in ode and check that y only shows as y(x) and not as just y #as the PDEtools:-Library:-GetDepVars([_self:-y],ode) code above does not detect this. #i.e. it does not check y'(x)+y=0 if numelems(indets(ode,identical(y))) > 0 then error("Parsing error, Can not have ",y," with no argument inside ",ode); fi; return y,x,the_order; end proc:
To use do
y,x,the_order := parse_single_ode(diff(y(x),x)+y(x)=0);
An alternative to the above is to pass the dependent function itself as well as the ode. This is what I do myself in my ode solver. Like this
parse_single_ode:=proc(ode::`=`,func::function(name))::symbol,symbol,integer; #parses single ode. returns back dep_var,indep_var,ode_order #throws exception when it detects parsing errors local y,x; local dep_variables_found::list,item; local the_order; if nops(func)<>1 then error("Parsing error, dependent variable must contain one argument, found ", func); fi; y:=op(0,func); x:=op(1,func); #basic verification if not has(ode,y) then error ("Supplied ode ",ode," has no ",y); fi; if not has(ode,x) then error ("Supplied ode ",ode," has no ",x); fi; if not has(ode,func) then error ("Supplied ode ",ode," has no ",func); fi; the_order := PDEtools:-difforder(ode,x); #note that the following call will also return y(x) if the input is not an ode #this will check that the dependent variable will show with SAME argument in the ode #i.e. if y(x) and y(t) show up in same ode, it will throw exception, which is what #we want. try dep_variables_found := PDEtools:-Library:-GetDepVars([y],ode); #print("dep_variables_found=",dep_variables_found); catch: error lastexception; end try; #now go over dep_variables_found and check the independent variable is same as x #i.e. ode can be y'(z)+y(z)=0 but function is y(x). for item in dep_variables_found do if not type(item,function) then error("Parsing error. Expected ",func," found ",item," in ode"); else if op(1,item) <> x then error("Parsing error. Expected ",func," found ",item," in ode"); fi; fi; od; #now go over all indents in ode and check that y only shows as y(x) and not as just y #as the PDEtools:-Library:-GetDepVars([_self:-y],ode) code above does not detect this. #i.e. it does not check y'(x)+y=0 if numelems(indets(ode,identical(y))) > 0 then error("Parsing error, Can not have ",y," with no argument inside ",ode); fi; return y,x,the_order; end proc:
To use do the same as before, but need to add \(y(x)\) as second argument. Like this
y,x,the_order := parse_single_ode(diff(y(x),x)+y(x)=0, y(x));
This is very much the same as above, but this function returns True if the ode is syntactly valid, else false. Can be used just to check if the ode is valid before using it.
It takes in the ode and the function \(y(x)\) and returns either true or false.
#March 15, 2024 export is_valid_single_ode :=proc(ode_in::`=`,func::function(name))::truefalse; local x:=op(1,func),y:=op(0,func); local ode:=ode_in; local dep_variables_found::list; local item; if nops(func)<>1 then RETURN(false); fi; if not has(ode,diff) then ode:=convert(ode,diff); fi; if `or`(not has(ode,diff), not has(ode,x), not has(ode,func)) then RETURN(false); fi; try dep_variables_found := PDEtools:-Library:-GetDepVars([y],ode); catch: RETURN(false); end try; map( X-> `if`( `or`(not type(X,function), op(1,X) <> x) ,RETURN(false),NULL), dep_variables_found ); #check there is no y on its own. Should always be y(x) if nops(indets(ode,identical(y))) <> 0 then RETURN(false); fi; RETURN(true); end proc:
Example usages
is_valid_single_ode(diff(y(x),x)+sin(x)=0, y(x) ) true
To test, do
L:=[[3*(D@@2)(y)(x)+diff(y(x),x)+1=sin(x),y(x),true], [diff(y(x),x)^2=z,y(x),true], [diff(y(x),x)^2=z,y(x),true], [diff(y(x),x)^2=y(x)^2,y(x),true], [diff(y(x),x)=y(z)^2,y(x),false], [1/diff(y(x),x)=y(x)^2,y(x),true], [y(x)=1,y(x),false], [diff(y(x),z)=1,y(x),false], [D(y)(x)=1,y(x),true], [(D@@2)(y)(x)+diff(y(x),x)=sin(x),y(x),true], [diff(y(x),x)+diff(y(z),z)=1,y(x),false], [diff(y(x),x)+diff(g(z),z)=1,y(x),true] ]: map(X->evalb(is_valid_single_ode(X[1],X[2])=X[3]),L); if andmap(x->evalb(x=true),%) then print("all tests passed"); else print("WARNING, not all tests passed"); fi; #gives [true, true, true, true, true, true, true, true, true, true, true, true] "all tests passed"
Suppose to want to check if an ode has \(y'(x)\) in it. We can not write has(ode, diff(y(x),x))
and see if this gives true or not, because this will also match \(y''\) and \(y'''\) and any higher
order.
One way is to first convert the ode to D form and then use has on D(y)(x). This will not
match \(y''\) anymore which is what we wanted, because \(y''\) becomes (D@@2)(y)(x) and so the check
for first order diff will work as expected.
ode:=diff(y(x),x$2)+3*y(x)=0; ode:=convert(ode,D); has(ode,(D)(y)(x)) #false
Another example
ode:=diff(y(x),x$2)+1/diff(y(x),x)=0; ode:=convert(ode,D); has(ode,(D)(y)(x)) #true
The same thing if we wanted to check if the ode has \(y''\) or not.
ode:=diff(y(x),x$3)+diff(y(x),x)=0; ode:=convert(ode,D); has(ode,(D@@2)(y)(x)) #false
ode:=diff(y(x),x$3)+ x*diff(y(x),x)+sin(x)=0; if has(DEtools:-odeadvisor(ode,y(x),['linear']),_linear) then print("linear"); else print("not linear ode"); fi
The above prints "linear"
ode:=y(x)*diff(y(x),x$3)+ x*diff(y(x),x)+sin(x)=0; if has(DEtools:-odeadvisor(ode,y(x),['linear']),_linear) then print("linear"); else print("not linear ode"); fi
The above prints "not linear ode". This works for any ode order.
ode:=diff(y(x),x$3)+ x*diff(y(x),x)+sin(x)=0; PDEtools:-difforder(ode)
Gives \(3\).
Given a linear ode of the form \(A y''+B y' + C y= f(x)\) how to find \(A,B,C,f(x)\) ?
ode:=diff(y(x),x$3)+ x*diff(y(x),x)+99*y(x)=sin(x); L:=DEtools:-convertAlg(ode,y(x)); #this only works on linear ode's
Gives L := [[99, x, 0, 1], sin(x)].
\(L\) is a list. The second entry is \(f(x)\) and the first entry of \(L\) is a list which gives the coefficents of the ode. Notice they are ordered from lowest order to highest order of the ode. Since this is third order ode, then there are \(4\) entries. The first is the coefficient of \(y(x)\), the second is the coefficient of \(y'\) and the third is the coefficient of \(y''\) and the fourth entry is the coefficient of \(y'''\). Notice that coefficient of \(y''\) is zero since \(y''\) is missing from the ode.
This function DEtools:-convertAlg only works on linear ode’s. Therefore we need to
check if the ode is linear first before using. How to check for linear ode is given
above.
I need to find the degree of the highest derivative in an ode.
For example, if the input is \[ x \left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right )^{5}+x \left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right )^{2}+\left (\frac {d}{d x}y \left (x \right )\right ) y \left (x \right )+\left (\frac {d^{3}}{d x^{3}}y \left (x \right )\right )^{4} \sin \left (x \right )+5 y \left (x \right )+\sin \left (x \right )+\frac {d^{6}}{d x^{6}}r \left (x \right )+\left (\frac {d^{3}}{d x^{3}}y \left (x \right )\right )^{4} \cos \left (x \right ) = 0 \]
Then the degree for highest derivative of \(y\) w.r.t. \(x\) is 4. For \[ y''(x)^2 + y'(x)+y(x)=0 \] It will be 2.
This function returns the order and degree of such term.
find_all_derivatives_of_specific_order:=proc(expr,y::symbol,x::symbol,N::posint)::set; local t1,t2; if not has(expr,diff(y(x),x$N)) then return {}; fi; t1 := identical(diff(y(x),x$N))^anything; t2 := identical(diff(y(x),x$N)); return indets[flat](expr,{t1,t2}); #MUST use flat end proc: get_order_and_degree_of_largest_derivative:=proc(expr,y::symbol,x::symbol)::integer,anything; local the_order,the_degree; local cand::set; local the_exponent; local item; if not has(expr,diff(y(x),x)) then return 0,0; fi; the_order := PDEtools:-difforder(expr,x); cand := find_all_derivatives_of_specific_order(expr,y,x,the_order); if nops(cand)=0 then the_degree := 0; else the_degree:=1; for item in cand do if type(item,`^`) then the_exponent := op(2,item); if type( the_exponent,symbol) or not type( the_exponent,numeric) then #assume larger the_degree :=the_exponent; else if type(the_degree,symbol) or not type( the_exponent,numeric) then next; else if op(2,item)>the_degree then the_degree := op(2,item); fi; fi; fi; else if type(the_degree,symbol) then next; else if the_degree=0 then the_degree := 1; fi; fi; fi; od; fi; return the_order,the_degree; end proc:
And now it can be called as follows
ode:=diff(y(x),x$2)+diff(y(x),x)+y(x)=0; get_order_and_degree_of_largest_derivative(lhs(ode),y,x) #2,1 ode:=diff(y(x),x$2)*diff(y(x),x$2)^4+diff(y(x),x)+y(x)=0; get_order_and_degree_of_largest_derivative(lhs(ode),y,x) #2,5 ode:=y(x)=0; get_order_and_degree_of_largest_derivative(lhs(ode),y,x) #0,0 ode:=diff(y(x),x)^n=0; get_order_and_degree_of_largest_derivative(lhs(ode),y,x) #1,n
Given \[ a \left (\frac {d}{d x}y \left (x \right )\right )+b \left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right )+x +\cos \left (x \right )+y \left (x \right )+c \left (\frac {d}{d x}y \left (x \right )\right )^{2} = \sin \left (x \right ) \] How to move all terms with derivative to LHS side and everything else to RHS?
ode:= a*diff(y(x),x)+b*diff(y(x),x$2)+x+cos(x)+y(x)+c*diff(y(x),x)^2=sin(x); ode:=lhs(ode)-rhs(ode); LHS,RHS:=selectremove(has,ode,'diff'); new_ode:=LHS=-RHS;
\[ a \left (\frac {d}{d x}y \left (x \right )\right )+b \left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right )+c \left (\frac {d}{d x}y \left (x \right )\right )^{2} = -x -\cos \left (x \right )-y \left (x \right )+\sin \left (x \right ) \]
I had a need to find all derivatives of form diff(y( anything), anything) in an ode so to
check that \(y\) argument is not different among them.
For an example, given \[ a \left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right ) \left (\frac {d^{3}}{d x^{3}}y \left (x \right )\right )-\sqrt {1+\left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right )^{2}}+\frac {1}{\sin \left (\frac {d^{5}}{d x^{5}}y \left (x \right )\right )} = \frac {d}{d z}y \left (z \right )+{\mathrm e}^{y \left (x \right )+\frac {d}{d x}y \left (x \right )}+\frac {d}{d x}r \left (x \right ) \]
the result should be \[ \left [\begin {array}{c} \frac {d^{5}}{d x^{5}}y \left (x \right ) \\ \frac {d^{4}}{d x^{4}}y \left (x \right ) \\ \frac {d^{3}}{d x^{3}}y \left (x \right ) \\ \frac {d^{2}}{d x^{2}}y \left (x \right ) \\ \frac {d}{d x}y \left (x \right ) \\ \frac {d}{d z}y \left (z \right ) \end {array}\right ] \]
One issuse is how to check for diff and also check that the dependent variable is \(y\) so as not to
pick other dependent variables such as \(z\) in this example. This was done by converting diff to
D otherwise it will not work.
restart; expr:=a*diff(y(x),x$2)*diff(y(x),x$3)-sqrt(1+ diff(y(x),x$2)^2)+1/sin(diff(y(x),x$5))=diff(y(z),z)+exp(y(x)+diff(y(x),x))+diff(r(x),x); #this finds all derivatives list_of_diffs:=indets(expr,'satisfies'(s->op(0,s)='diff' and op([0,1],convert(s,D))=y)); #This finds all dependent variables list_of_diffs:=convert(list_of_diffs,list); list_of_indep_variables := map(x->PDEtools:-Library:-GetIndepVars(x)[-1],list_of_diffs); #this converts it to set. If the ODE is valid, then the list_of_indep_variables should #have one entry $x$ in it and nothing else. list_of_indep_variables := convert(ListTools:-Flatten(list_of_indep_variables),set); if nops(list_of_indep_variables)>1 then error( cat("Only one independent variable expected in differential form, found ", convert(list_of_indep_variables,string)) ); fi; if list_of_indep_variables[1]<>x then error( cat("Independent variable expected in differential form not same as independent variable in function given ", convert(list_of_indep_variables,string)) ); fi;
Another option instead of doing all the above is to do this
expr:=a*diff(y(x),x$2)*diff(y(x),x$3)-sqrt(1+ diff(y(x),x$2)^2)+1/sin(diff(y(x),x$5))=diff(y(z),z)+exp(y(x)+diff(y(x),x))+diff(r(x),x); try PDEtools:-Library:-GetDepVars([y(x)],expr); catch error "functions with name [y] having different dependency: [[y(x), y(z)]]" end try;
The function PDEtools:-Library:-GetDepVars([y(x)],expr) checks that only \(y(x)\)
dependency shows up. It throws an error otherwise. So if an error is thrown, then this means
\(y\) shows up with different independent variables.
Sometimes it is useful to invert an ode. i.e. make the independent variable the dependent variable, and the dependent variable the independent. For example, given \[ 1+\left (\frac {x}{y \left (x \right )}-\sin \left (y \left (x \right )\right )\right ) \left (\frac {d}{d x}y \left (x \right )\right ) = 0 \] We want the ode to become \[ -\sin \left (y \right ) y +y \left (\frac {d}{d y}x \left (y \right )\right )+x \left (y \right ) = 0 \]
This can be done as follows
restart; ode:=1+ (x/y(x)-sin(y(x) ))*diff(y(x),x)=0; tr:={x=u(t),y(x)=t}; ode:=PDEtools:-dchange(tr,ode); ode:=eval(ode,[t=y,u=x]); ode:=simplify(ode);
\[ \frac {-\sin \left (y \right ) y +y \left (\frac {d}{d y}x \left (y \right )\right )+x \left (y \right )}{y \left (\frac {d}{d y}x \left (y \right )\right )} = 0 \]
In this case, we can get rid of the denomator, but this is a manual step for now.
ode:=numer(lhs(ode))=0;
\[ -\sin \left (y \right ) y +y \left (\frac {d}{d y}x \left (y \right )\right )+x \left (y \right ) = 0 \]
The above can now be solved more easily for \(x(y)\) than solving the orignal ode for \(y(x)\).
For say Bessel ode of order zero:
eq:= x^2*diff(y(x),x$2)+x*diff(y(x),x)+x^2*y(x)=0; DEtools[indicialeq](eq,x,0,y(x)); #x^2 = 0
The third argument above is the singularity point of interest. So we have two roots, both zero. These are now used for finding the power series solution \(y(x)\) if needed.
Another example, is Bessel of order 1
eq:= x^2*diff(y(x),x$2)+x*diff(y(x),x)+(x^2-1)*y(x)=0; DEtools[indicialeq](eq,x,0,y(x)); #x^2-1 = 0
To write \(y'(x)=x\), one way is diff(y(x),x)=x and another is D(y)(x)=x. To write \(y''(x)=x\), one way is
diff(y(x),x$2)=x and another is (D@@2)(y)(x)=x.
To convert from one form to another use convert(eq,diff) or convert(eq,D)
to solve \(\frac {\partial u(x,t)}{\partial t}=k \frac {\partial ^2 u(x,t)}{\partial x^2}\) with homogeneous dirichlet boundary conditions \(u(0,t)=0,u(L,t)=0\) the commands are
restart; pde:=diff(u(x,t),t)=k*diff(u(x,t),x$2); bc:=u(0,t)=0,u(L,t)=0; sol:=pdsolve([pde,bc]) assuming 0<L:
Which gives
\[ u \left ( x,t \right ) =\sum _{{\it \_Z1}=1}^{\infty }{\it \_C1} \left ( {\it \_Z1} \right ) \sin \left ( {\frac {\pi \,{\it \_Z1}\,x}{L}} \right ) {{\rm e}^{-{\frac {k{\pi }^{2}{{\it \_Z1}}^{2}t}{{L}^{2}}}}} \]
Which can be made more readable as follows
sol:=algsubs(_Z1=n,sol): sol:=algsubs(Pi*n/L=lambda(n),sol);
\[ u \left ( x,t \right ) =\sum _{n=1}^{\infty }{\it \_C1} \left ( n \right ) \sin \left ( x\lambda \left ( n \right ) \right ) {{\rm e}^{-kt \left ( \lambda \left ( n \right ) \right ) ^{2}}} \]
For homogeneous Neumann B.C., at \(x=0\), let \(\frac {\partial u}{\partial x}=0\) and at \(x=L\) let \(u(L,t)=0\), the solution it gives looks different than my hand solution
restart; pde:=diff(u(x,t),t)=k*diff(u(x,t),x$2); bc:=D[1](u)(0,t)=0,u(L,t)=0; pdsolve([pde,bc]) assuming 0<L;
It gives
\[ u \left ( x,t \right ) ={\it \_C3}\,{\it \_C2}\, \left ( {{\rm e}^{1/4\,{\frac {2\,i\pi \,xL-k{\pi }^{2}t}{{L}^{2}}}}}+{{\rm e}^{-1/4\,{\frac {\pi \, \left ( 2\,ixL+k\pi \,t \right ) }{{L}^{2}}}}} \right ) \]
I need to look more into the above and see if this comes out to be the same as my hand solution.
Another example, with initial conditions now given
restart; pde:=diff(u(x,t),t)=k*diff(u(x,t),x$2); bc:=D[1](u)(0,t)=0,u(L,t)=0; ic:=u(x,0)=f(x); sol:=pdsolve([pde,bc,ic],u(x,t)) assuming 0<L; sol1:=algsubs(_Z2=n,sol);
The result is
\[ u \left ( x,t \right ) =\sum _{n=1}^{\infty } \left ( 2\,{\frac {1}{L}{{\rm e}^{-1/4\,{\frac {k{\pi }^{2}t \left ( 1+2\,n \right ) ^{2}}{{L}^{2}}}}}\cos \left ( 1/2\,{\frac {\pi \,x \left ( 1+2\,n \right ) }{L}} \right ) \int _{0}^{L}f \left ( x \right ) \cos \left ( 1/2\,{\frac {\pi \,x \left ( 1+2\,n \right ) }{L}} \right ) \,{\rm d}x} \right ) \]
Another example
restart; pde:=diff(u(x,t),t)=k*diff(u(x,t),x$2); bc:=D[1](u)(0,t)=0,u(L,t)=0; ic:=u(x,0)=3*sin(Pi*x/L)-sin(3*Pi*x/L); sol:=pdsolve([pde,bc,ic],u(x,t)) assuming 0<L; sol1:=algsubs(_Z2=n,sol);
\[ u \left ( x,t \right ) =\sum _{n=1}^{\infty }768\,{\frac {1}{\pi \, \left ( 16\,{n}^{4}+32\,{n}^{3}-136\,{n}^{2}-152\,n+105 \right ) }{{\rm e}^{-1/4\,{\frac {k{\pi }^{2}t \left ( 1+2\,n \right ) ^{2}}{{L}^{2}}}}}\cos \left ( 1/2\,{\frac {\pi \,x \left ( 1+2\,n \right ) }{L}} \right ) } \]
Another example
restart; pde:=diff(u(x,t),t)=k*diff(u(x,t),x$2); bc:=u(0,t)=0,u(L,t)=0; ic:=u(x,0)=3*sin(Pi*x/L)-sin(3*Pi*x/L); sol:=pdsolve([pde,bc,ic],u(x,t)) assuming 0<L;
\[ u \left ( x,t \right ) =\sin \left ( {\frac {\pi \,x}{L}} \right ) {{\rm e}^{-9\,{\frac {{\pi }^{2}kt}{{L}^{2}}}}} \left ( -2\,\cos \left ( 2\,{\frac {\pi \,x}{L}} \right ) +3\,{{\rm e}^{8\,{\frac {{\pi }^{2}kt}{{L}^{2}}}}}-1 \right ) \]
The above answer seems wrong. There is not even a summation in it. It is different from my hand solution. Look more into it.
Add this
expr:=diff(y(x),x); Typesetting:-Settings(typesetprime=true, prime=x):
The above will display the expression as \(y'(x)\). To make it now show the \(x\) do
expr:=diff(y(x),x); Typesetting:-Settings(typesetprime=true, prime=x): Typesetting:-Suppress(y(x));
Now it will show the expression as just \(y'\). For all the above to work, make sure you have
Typesetting level set to Extended in the GUI.
This is done inside Tools->Options->Display menu.
To clear all the above Typesetting, do restart or do Typesetting:-Unsuppress(y(x))
The Maple syntax for seeting initial and boundary conditions is very confusing, as compared to Mathematica, which seems to me to be simpler. So I wrote this to remind me of the syntax each time.
For PDE, assuming dependent variable is \(u(x,t)\) then
| Conditions | Maple code |
| \(u(0,t)=0\) | u(0,t)=0 |
| \(\frac {\partial u}{\partial x}=0\) at \(x=0\) | D[1](u)(0,t)=0 |
| \(\frac {\partial ^2 u}{\partial x^2}=0\) at \(x=0\) | D[1,1](u)(0,t)=0 |
| \(\frac {\partial ^3 u}{\partial x^3}=0\) at \(x=0\) | D[1,1,1](u)(0,t)=0 |
| \(\frac {\partial u}{\partial t}=0\) at \(t=0\) | D[2](u)(x,0)=0 |
| \(\frac {\partial ^2 u}{\partial t^2}=0\) at \(t=0\) | D[2,2](u)(x,0)=0 |
| \(\frac {\partial ^3 u}{\partial t^3}=0\) at \(t=0\) | D[2,2,2](u)(x,0)=0 |
Notice the syntax for the last one above. It is (D[1]@@2)(u)(0,t)=0 and not
(D@@2)[1](u)(0,t)=0
For an ODE, assuming dependent variable is \(y(x)\) then the syntax is
| Conditions | Maple code |
| \(y(0)=0\) | y(0)=0 |
| \(\frac {dy}{dx}=0\) at \(x=0\) | D(y)(0)=0 |
| \(\frac {d^2 y}{d x^2}=0\) at \(x=0\) | (D@@2)(y)(0)=0 |
Use dsolve(ode,Lie)
To find symmetries, do
DEtools:-symgen(ode,y(x),HINT=[c__1+c__2*x+c__3*y,c__4+c__5*x+c__6*y])
or just
DEtools:-symgen(ode,y(x))
To debug it do
stopat(`ODEtools/symgen`); before calling dsolve or DEtools:-symgen
given an ode
\[ y \left (x \right ) = \left (\frac {d}{d x}y \left (x \right )\right )^{3} y \left (x \right )^{2}+2 x \left (\frac {d}{d x}y \left (x \right )\right ) \]
do change of variable \(u(x)=y(x)^2\)
restart; ode:=y(x)=diff(y(x),x)^3*y(x)^2+2*x*diff(y(x),x); new_ode:=PDEtools:-dchange({y(x)=sqrt(u(x))},ode,{u});
\[ \sqrt {u \left (x \right )} = \frac {\left (\frac {d}{d x}u \left (x \right )\right )^{3}}{8 \sqrt {u \left (x \right )}}+\frac {x \left (\frac {d}{d x}u \left (x \right )\right )}{\sqrt {u \left (x \right )}} \]
given an ode
\[ \frac {d^{2}}{d t^{2}}y \left (t \right )+y \left (t \right ) = 2 t \]
do change of variable \(t=\tau + \pi \)
restart; ode:=diff(y(t),t$2)+y(t)=2*t; PDEtools:-dchange({t=tau+Pi},ode,known={t},unknown={tau},params=Pi)
\[ \frac {d^{2}}{d \tau ^{2}}y \left (\tau \right )+y \left (\tau \right ) = 2 \tau +2 \pi \]
it is important to use params=Pi. Watch what happens if we do not do that
restart; ode:=diff(y(t),t$2)+y(t)=2*t; PDEtools:-dchange({t=tau+Pi},ode,known={t},unknown={tau});
\[ \frac {d^{2}}{d \tau ^{2}}y \left (\tau , \pi \right )+y \left (\tau , \pi \right ) = 2 \tau +2 \pi \] Which is not what we want.
This verifies solution given in https://math.stackexchange.com/questions/3477732/can-t-see-that-an-ode-is-equivalent-to-a-bessel-equation Where a change of variables on \[ \xi ^2 \frac {d^2\eta }{d\xi ^2} + \xi \frac {d\eta }{d\xi } - (\xi ^2+n^2)\eta =0 \] Was made using \[ \eta =\frac {y}{x^\alpha }, \quad \xi =\beta x^\gamma , \] To produce the ode \[ \frac {d^2y}{dx^2} - \frac {(2\alpha -1)}{x}\frac {dy}{dx} - (\beta ^2 \gamma ^2 x^{2\gamma -2}+\frac {n^2\gamma ^2-\alpha ^2}{x^2})y=0. \] In Maple the above is done using
restart; ode := zeta^2*diff(eta(zeta),zeta$2) + zeta*diff(eta(zeta),zeta) - (zeta^2 + n^2)*eta(zeta) = 0; the_tr:={zeta=beta*x^gamma,eta(zeta)=y(x)/x^alpha}; PDEtools:-dchange(the_tr,ode,{y(x),x},'known'={eta(zeta)},'uknown'={y(x)},'params'={alpha,beta,gamma}); simplify(%); numer(lhs(%))=0; simplify(numer(lhs(%))/(x^(1-alpha)))=0; numer(lhs(%))=0; collect(%,[y(x),diff(y(x),x),diff(y(x),x$2)]);
Which gives \[ \left (-\gamma ^{2} x^{-1+2 \gamma } \beta ^{2} x -\gamma ^{2} n^{2}+\alpha ^{2}\right ) y \left (x \right )+\left (-2 \alpha x +x \right ) \left (\frac {d}{d x}y \left (x \right )\right )+x^{2} \left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right ) = 0 \]
Here is another example. Here to make change of variables to polar coordinates by making \(x=r \cos \theta \) and \(y=r \sin \theta \) The ode is \[ \frac {y-x y'}{\sqrt {1+(y')^2}}=x^2+y^2 \]
In Maple
restart; ode := (y(x)-x*diff(y(x),x))/sqrt(1+ diff(y(x),x)^2) = x^2+y(x)^2; the_tr:={x=r(t)*cos(t),y(x)=r(t)*sin(t)}; PDEtools:-dchange(the_tr,ode,{r(t),t},'known'={y(x)},'uknown'={r(t)});
Which gives
\[ \frac {r \left (t \right ) \sin \left (t \right )-\frac {r \left (t \right ) \cos \left (t \right ) \left (\left (\frac {d}{d t}r \left (t \right )\right ) \sin \left (t \right )+r \left (t \right ) \cos \left (t \right )\right )}{\left (\frac {d}{d t}r \left (t \right )\right ) \cos \left (t \right )-r \left (t \right ) \sin \left (t \right )}}{\sqrt {1+\frac {\left (\left (\frac {d}{d t}r \left (t \right )\right ) \sin \left (t \right )+r \left (t \right ) \cos \left (t \right )\right )^{2}}{\left (\left (\frac {d}{d t}r \left (t \right )\right ) \cos \left (t \right )-r \left (t \right ) \sin \left (t \right )\right )^{2}}}} = r \left (t \right )^{2} \left (\cos ^{2}\left (t \right )\right )+r \left (t \right )^{2} \left (\sin ^{2}\left (t \right )\right ) \]
Here is another example. Where we want to change \(R(r)\) to \(y(x)\) everywhere
\[ \frac {d^{2}}{d r^{2}}R \left (r \right )+\frac {d}{d r}R \left (r \right )+R \left (r \right ) = 0 \]
restart; ode:=diff(R(r),r$2)+diff(R(r),r)+R(r)=0; the_tr:={r=x,R(r)=y(x)}; PDEtools:-dchange(the_tr,ode,{y(x),x},'known'={R(r),r},'uknown'={y(x),x});
\[ \frac {d^{2}}{d x^{2}}y \left (x \right )+\frac {d}{d x}y \left (x \right )+y \! \left (x \right ) = 0 \]
The format of the transformation is old_independent_variable=new_independent_variable
and old_dependent_variable=new_dependent_variable
Make a phase plot of \[ y'(x) = \sqrt {y(x)^{2}-1} \]
The phase plot has \(x\) on the x axis and has \(y\) on the y axis. It shows the family of solutions for different initial conditions.
restart; ode:=diff(y(x), x) = sqrt(y(x)^2 - 1); DEtools:-DEplot(ode,y(x),x=-2..2,y=1..3)
To show specific solution curve that passes via some initial conditions such as \(y(0)=2\) then do
restart; ode:=diff(y(x), x) = sqrt(y(x)^2 - 1); DEtools:-DEplot(ode,y(x),x=-2..2,y=1..3,[y(0)=2])
Make a phase plot of \[ \frac {d^{2}}{d t^{2}}x \left (t \right )+\frac {\frac {d}{d t}x \left (t \right )}{2}+x \left (t \right ) = u \left (t \right ) \]
By plotting \(x(t)\) vs \(x'(t)\) without solving the ODE.
restart; alias(DS=DynamicSystems): ode := diff(x(t),t$2) +1/2*diff(x(t),t)+ x(t) = u(t); sys:=DS:-DiffEquation(ode,'outputvariable'=[x(t)],'inputvariable'=[u(t)]); sys0:=DS:-StateSpace(sys); eq1:=diff(x1(t),t)=sys0:-a[1,..].Vector([x1(t),x2(t)]); eq2:=diff(x2(t),t)=sys0:-a[2,..].Vector([x1(t),x2(t)]); DEtools:-DEplot([eq1,eq2],[x1(t),x2(t)],t=0..35,[[x1(0)=1,x2(0)=1]],x1=-2..2,x2=-2..2, numpoints=200, linecolor=black, axes=boxed);
Given ode where \(y\) is dependent variable and \(x\) is independent variable, how to normalize the ode so that all terms with \(y(x)\) are on left side and everything is on the RHS? THis makes it easier to see what the forcing function is. i.e. we want to write the ode as
\[ a y''+b y'+c y= f(x) \]
Here is function to do this. This works for any ode or any expression.
export move_y_x_to_each_side:=proc(expr_in::`=`,y::symbol,x::symbol)::`=`; local expr:=expr_in; local stuff_with_y,stuff_not_with_y; expr := lhs(expr)-rhs(expr); if not expr::`+` then RETURN(expr_in); fi; expr := numer(normal(expr)); stuff_with_y,stuff_not_with_y := selectremove(has,expr,y); stuff_not_with_y := convert(stuff_not_with_y,diff); RETURN(stuff_with_y = -stuff_not_with_y); end proc:
Examples
ode:=diff(y(x),x$2)+x*y(x)+sin(x)-3=tan(x); expr:=move_y_x_to_each_side(ode,y,x); # diff(y(x), x, x) + x*y(x) = 3 - sin(x) + tan(x)
And
ode:=diff(y(x),x$2)+x*y(x)+sin(x)-3=tan(x)+y(x)+Pi; expr:=move_y_x_to_each_side(ode,y,x); # diff(y(x), x, x) + x*y(x) - y(x) = 3 - sin(x) + tan(x) + Pi
And
ode:=move_y_x_to_each_side(diff(x(t),t$2)/x(t) + diff(x(t),t)/x(t) - F(t)/x(t) =1,x,t) expr:=move_y_x_to_each_side(ode,x,t); # diff(x(t), t, t) + diff(x(t), t) - x(t) = F(t)
Sometimes I make solution for an ode with constants of integrations. These must be of form
_Cn where \(n\) is positive integer, such as _C1 or _C5. Now to integrate this solution
once more, need to come up with new constant that does not already show in the
solution.
To find current constants in an expression use this
sol:=x-_C1*x+_C3*x^2; indets(sol,And(symbol, suffixed(_C, nonnegint))); {_C1, _C3}
Therefore, one way to make a new constant of integration, is to find the largest numbered one and increase by one. Like this
restart; sol:=x-_C1*x+_C3*x^2; myconstants:=indets(sol,And(symbol, suffixed(_C, nonnegint))); map(X->String(X),myconstants); map(X->X[3..],%); map(X->:-parse(X),%); n:=max(%); new_constant:=_C||(n+1);
Sometimes it is usefull to know if ode is missing \(y(x)\) as this allows one to do the substitution \(y'=u\) and reduce the order of the ode.
This is a function which takes in an ode and returns true if it is missing \(y(x)\) or false otherwise.
#added 2/18/2024. checks if an ode is missing y export is_ode_missing_y:=proc(ode_in::`=`,y::symbol,x::symbol)::truefalse; local ode:=ode_in; local ode_order::posint := PDEtools:-difforder(ode_in); local N::posint; ode:=lhs(ode)-rhs(ode); ode:=convert(ode,D); for N from 1 to ode_order do ode:=eval(ode,[(D@@N)(y)(x)=Y||N]); od; RETURN(not has(ode,y(x))); end proc;
To use do
restart; ode:=diff(y(x),x$2)+diff(y(x),x$4)+x=0; is_ode_missing_y(ode,y,x); #true ode:=diff(y(x),x$2)+diff(y(x),x$4)+x=sin(x)+y(x); is_ode_missing_y(ode,y,x); #false
See http://www.maplesoft.com/applications/view.aspx?SID=1533&view=html&L=G
use select. For example
>restart; >my_list:=[1,3.4,3+I,5]; >select(x->evalb(Im(x)=0),my_list); [1, 3.4, 5]
restart; m:=Matrix( [[1.3,2,3],[3,4,4] ]); matrixTestQ := proc(m::Matrix) local r,c,i,j; (r,c):=LinearAlgebra[Dimensions](m); for i from 1 to r do for j from 1 to c do if( not evalb( whattype(m[i,j]) = integer) ) then return(false); end if; end do; end do; return true; end proc; >matrixTestQ(m); false
I am sure there is a better way than the above. Need to find out.
Given
\[ 3+x +\sqrt {-4 a c +b^{2}}+\sin \left (y \right )+x^{3} \sqrt {39}+\sqrt {\cos }x \] Find terms that are sqrt. Use indets
restart; expr_with_radical:= 3+x+sqrt(b^2-4*a*c)+sin(y)+x^3*sqrt(39)+sqrt(cos(x));| indets(expr_with_radical, algebraic^fraction)
\( \{ \sqrt {39}, \sqrt {-4 a c +b^{2}}, \sqrt {\cos }x \} \)
Alternative is to use type radical
restart; expr_with_radical:= 3+x+sqrt(b^2-4*a*c)+sin(y)+x^3*sqrt(39)+sqrt(cos(x));| indets(expr_with_radical, radical)
\( \{ \sqrt {39}, \sqrt {-4 a c +b^{2}}, \sqrt {\cos }x \} \)
I wanted to simplify an expression which could have csgn() in it, and find all the
arguments.
\[ \frac {1+\mathrm {csgn} \left (a \right ) a}{3 \mathrm {csgn} \left (b \right ) b} \] One way is
restart; expr:=(1+csgn(a)*a)/(3*csgn(b)*b): expr:=subsindets(expr,'specfunc( anything, csgn )',f->1);
\[ \frac {1+a}{3 b} \]
Given sol:=1/2*2^(1/2)*csgn(x)*x*csgn(y); how to find all symbols inside csgn which
will be \(x,y\) in this case?
restart; sol:=1/2*2^(1/2)*csgn(x)*x*csgn(y); indets(sol,'specfunc( anything, csgn )'); vars:=subsindets(%,'specfunc( anything, csgn )',f->op(f))
Gives {x, y}
Now if we want to simplify the above solution by assuming that all variables inside vars are
positive, how to do that?
restart; sol:=1/2*2^(1/2)*csgn(x)*x*csgn(y); indets(sol,'specfunc( anything, csgn )'); vars:=subsindets(%,'specfunc( anything, csgn )',f->op(f)); simplify(sol) assuming op(map2(`<`,0,vars))
Gives \(\frac {\sqrt {2}\, x}{2}\). Notice in the above the use of op(map2(`<`,0,vars)), this will generate the
sequence 0 < x, 0 < y automatically. op is needed otherwise the result will be
{0 < x, 0 < y} which will give syntax error when passed to assuming
Ofcourse, it would have been also possible to just write
simplify(sol) assuming positive;
And get the same result. But sometimes we might want to specify which variables are to be assumed positive and not all of them at once in the expression.
I wanted to replace |expr| by (expr)
One way is
restart; expr:=u(x) = _C1*exp(-3*x^(1/3)*sqrt(c))*(3*x^(1/3)*sqrt(c) + 1) + _C2*exp(3*x^(1/3)*sqrt(c))*abs(-1 + 3*x^(1/3)*sqrt(c));
\[ u \left (x \right ) = \textit {\_C1} \,{\mathrm e}^{-3 x^{\frac {1}{3}} \sqrt {c}} \left (3 x^{\frac {1}{3}} \sqrt {c}+1\right )+\textit {\_C2} \,{\mathrm e}^{3 x^{\frac {1}{3}} \sqrt {c}} {| -1+3 x^{\frac {1}{3}} \sqrt {c}|} \]
expr:=subsindets(expr,'specfunc( anything, abs )',f->op(f));
\[ u \left (x \right ) = \textit {\_C1} \,{\mathrm e}^{-3 x^{\frac {1}{3}} \sqrt {c}} \left (3 x^{\frac {1}{3}} \sqrt {c}+1\right )+\textit {\_C2} \,{\mathrm e}^{3 x^{\frac {1}{3}} \sqrt {c}} \left (-1+3 x^{\frac {1}{3}} \sqrt {c}\right ) \]
For an example, How to find list of all \(\ln \) functions in this expression?
\[ \ln \left ({| x +1|}\right )+2 x \ln \left (x \right )+\sin \left (x \right ) \]
restart; expr:=ln(abs(x+1))+2*x*ln(x)+sin(x); tmp := indets(expr,'specfunc(anything,ln)'); # tmp := {ln(x), ln(abs(x + 1))}
To pick only \(\ln \) functions which has \(abs\) inside them anywhere, replace the above with
restart; expr:=ln(abs(x+1))+2*x*ln(x)+sin(x); lis:=indets(expr,'specfunc(anything,ln)'); select(Z->has(Z,abs),lis) # tmp := {ln(abs(x + 1))}
Or, better alternative to the above is
restart; expr:=ln(abs(x+1))+2*x*ln(x)+sin(x); indets(expr,'specfunc( satisfies(u->has(u,abs)) ,ln )'); # tmp := {ln(abs(x + 1))}
Given
\[ \sin \left (x \right )+\ln \left ({| x |}\right )+\ln \left (x +\frac {{| y |}}{\sqrt {{| x +3|}}}\right )+\ln \left (x^{3}\right )+\cos \left ({| x |}\right ) \]
How to remove the absolute, the ones only inside each \(\ln \) in the above expression?
restart; expr:=sin(x)+ln(abs(x))+ln(x+abs(y)/sqrt(abs(x+3)))+ln(x^3)+cos(abs(x)); expr:=evalindets(expr,'specfunc(ln)',f->evalindets(f,'specfunc(abs)',f->op(1,f))) # sin(x) + ln(x) + ln(x + y/sqrt(x + 3)) + ln(x^3) + cos(abs(x))
\[ \sin \left (x \right )+\ln \left (x \right )+\ln \left (x +\frac {y}{\sqrt {x +3}}\right )+\ln \left (x^{3}\right )+\cos \left ({| x |}\right ) \]
Given
\[ -\frac {\left (\ln \left (\frac {\left (b +\sqrt {b^{2}+y \left (x \right )^{2}}\, \mathrm {signum}\left (b \right )\right ) b}{y \left (x \right )}\right )+\ln \left (2\right )\right ) \mathrm {signum}\left (b \right )}{b} = \mathit {\_C1} +\frac {-\ln \left (a \right )+\ln \left (x \right )-\ln \left (a +\sqrt {a^{2}+x^{2}}\, \mathrm {signum}\left (a \right )\right )-\ln \left (2\right )}{{| a |}} \]
How to find all arguments of signum and simplify the above by assuming they are all
positive?
restart; expr:=-(ln((b + sqrt(b^2 + y(x)^2)*signum(b))*b/y(x)) + ln(2))*signum(b)/b = _C1 + (-ln(a) + ln(x) - ln(a + sqrt(a^2 + x^2)*signum(a)) - ln(2))/abs(a); lis:=indets(expr,'specfunc(anything,signum)'); assum:=convert(map(x->op(1,x)>0,lis),list); simplify(expr,assume=assum);
\[ \frac {-\ln \left (b \right )-\ln \left (\frac {b +\sqrt {b^{2}+y \left (x \right )^{2}}}{y \left (x \right )}\right )-\ln \left (2\right )}{b} = \frac {\mathit {\_C1} a -\ln \left (a \right )-\ln \left (a +\sqrt {a^{2}+x^{2}}\right )+\ln \left (x \right )-\ln \left (2\right )}{a} \]
Given
\[ -\frac {\left (\ln \left (\frac {\left (b +\sqrt {b^{2}+y \left (x \right )^{2}}\, \mathrm {signum}\left (b \right )\right ) b}{y \left (x \right )}\right )+\ln \left (2\right )\right ) \mathrm {signum}\left (b \right )}{b} = \mathit {\_C1} +\frac {-\ln \left (a \right )+\ln \left (x \right )-\ln \left (a +\sqrt {a^{2}+x^{2}}\, \mathrm {signum}\left (a \right )\right )-\ln \left (2\right )}{{| a |}} \]
How to replace all signum by 1?
restart; expr:=-(ln((b + sqrt(b^2 + y(x)^2)*signum(b))*b/y(x)) + ln(2))*signum(b)/b = _C1 + (-ln(a) + ln(x) - ln(a + sqrt(a^2 + x^2)*signum(a)) - ln(2))/abs(a); evalindets(expr, 'specfunc(anything,signum)', f -> 1);
\[ -\frac {\ln \left (\frac {\left (b +\sqrt {b^{2}+y \left (x \right )^{2}}\right ) b}{y \left (x \right )}\right )+\ln \left (2\right )}{b} = \textit {\_C1} +\frac {-\ln \left (a \right )+\ln \left (x \right )-\ln \left (a +\sqrt {a^{2}+x^{2}}\right )-\ln \left (2\right )}{{| a |}} \]
Given expression \(3 \sin \left (x \right )+t +3 f \left (x , t\right ) t +g \left (x , t\right )\) find if it contains function \(f()\).
Use indets with specfunc(f)
restart; expr := 3*sin(x)+t+3*f(x,t)*t+g(x,t); res := indets(expr, specfunc(f)); if numelems(res)<>0 then print("Found f(x,t)"); else print("could not find f(x,t)"); fi;
"Found f(x,t)"
Given expression \(3 \sin \left (x \right )+t +3 f \left (x , t\right ) t +g \left (x , t\right )\) find all functions, if any, in the expression.
Use indets with function
restart; expr := 3*sin(x)+t+3*f(x,t)*t+g(x,t); res := indets(expr,function); if numelems(res)<>0 then print("Found these functions",res); else print("could not find any function)"); fi;
"Found these functions", f(x, t), g(x, t), sin(x)
Given expression \(3 \sin \left (x \right )+t +3 f \left (x , t\right ) t +g \left (x , t\right )\) find all functions, if any, in the expression but exclude the math functions such as \(\sin \) in the above.
restart; expr := 3*sin(x)+t+3*f(x,t)*t+g(x,t); res := indets(expr, And( function, Not(typefunc(mathfunc)))); if numelems(res)<>0 then print("Found these functions",res); else print("could not find any function)"); fi;
"Found these functions", f(x, t), g(x, t)
use op
restart; op(1..,f(x,t)) x, t
Note that op(0,f(x,t)) finds the function name.
restart; expr := 3*sin(z)+t+3*f(z,t,y)*t+g(x,t); res := indets(expr, patfunc(identical(z), anything)); if numelems(res)<>0 then print("Found these functions",res); else print("could not find any function)"); fi;
gives
"Found these functions", f(z, t, y), sin(z)
expr := 3*sin(z)+t+3*f(z,t,y)*t+g(x,t); res := indets(expr, patfunc(anything, identical(t), anything)); if numelems(res)<>0 then print("Found these functions",res); else print("could not find any function)"); fi;
gives
"Found these functions", f(z, t, y), g(x, t)
Given expression such as \(3+(1+x)\sin x\) or \(3+(1+x)\sin ^2 x\) use select to find any polynomial * sin^n subexpressions.
restart; mytype_1 := ''`*`'( {polynom(And(algebraic, satisfies(u -> not has(u, I))),x), Or( 'specfunc(sin)'^integer, 'specfunc(sin)') } )'; select(type, 3+(1+x)*sin(x),mytype_1); select(type, 3+(1+x)*sin(x)^2,mytype_1);
Gives
(1 + x) sin(x) 2 (1 + x) sin(x)
restart; my_type:=''`*`'( { Or('specfunc'(sin),'specfunc'(sin)^Or(integer,rational)), Or('specfunc'(cos),'specfunc'(cos)^Or(integer,rational))})'; type(sin(x)^2*cos(x)^3,my_type); type(sin(x)^2*cos(x),my_type); type(sin(x)*cos(x),my_type); type(cos(x)*sin(x)^(1/2),my_type); true true true true
I could not find a way to avoid writing Or('specfunc'(sin),'specfunc'(sin)^Or(integer,rational)
in order to match both \(\sin x\) and \(\sin ^2 x\). For these things, I find Mathematica patterns more flexiable.
The above can be done as follows in Mathematica
ClearAll[x,n,m,any] patt=any_.*Sin[_]^n_. * Cos[_]^m_. MatchQ[Sin[x]^2*Cos[2*x]^3,patt] MatchQ[Sin[x]^2*Cos[x],patt] MatchQ[Sin[x]*Cos[x],patt] MatchQ[Cos[x]*Sin[x],patt] True True True True
In Mathematica n_. says basically to match \(\sin x\) or \(\sin ^2 x\) since the dot says to match zero or
more. So no need to duplicate things as I did above in Maple. There might be a
way to do the same in Maple using structured type, but I could not find it. In
General, I find patterns in Mathematica more flexible and easier to use for this sort
of thing. Maple has patmatch command, but not as easy to use as Patterns in
Mathematica.
use indets with type 'indexed'
expr:=16*a[3]+6*a[1]; terms:=indets(expr,'indexed'); terms := {a[1], a[3]} #to find maximum index, then do map(x->op(x),terms) {1, 3}
Given say \(\frac {d^{2}}{d x^{2}}y \left (x \right )+n \left (\frac {d}{d x}y \left (x \right )\right )+3 = \sin \left (x \right )\) how to find all variables and functions in it, not including math functions such as \(\sin x\)?
So the result should be \(n,x,y(x)\).
ode:=diff(y(x),x$2)+n*diff(y(x),x)+3=sin(x); vars:=indets(ode, Or( And(symbol,Not(constant)), And(function,Not(typefunc(mathfunc)) ) )) #gives # vars := {n, x, diff(y(x), x), diff(y(x), x, x), y(x)}
I still need to work on excluding derivatives from the search.
I had case where I needed to check if something is integer or not. The problem is that the result had a symbol \(n\) in it. I need a way to tell Maple that to check if the result can be an integer given that \(n\) is also an integer.
Using type does not work, since can’t use assumptions. One way is to use coulditbe as
follows
restart; expr:=n-1+2*m; vars:=indets(expr,And(symbol,Not(constant))); coulditbe(expr,integer) assuming op(map(Z->Z::integer,vars)) # true
In the above indets(expr,And(symbol,Not(constant))) picks all variables in the
expression, and assuming op(map(Z->Z::integer,vars)) makes assumption that each is
integer.
Use
restart; expr:=4*Pi+sin(x); indets(expr,And(name,constant)) {Pi}
Given an expression such as
\[ \sin \left (x \right )+\left (\frac {d}{d x}y \left (x \right )\right )^{3}+\left (\frac {d}{d x}y \left (x \right )\right ) \ln \left (y \left (x \right ) \left (\frac {d}{d x}y \left (x \right )\right )^{2}\right )+\sqrt {\frac {d}{d x}y \left (x \right )}+\frac {x}{\left (\frac {d}{d x}y \left (x \right )\right )^{7}} \]
Find all \(y'(x)\) for any power that show up, so the result should be
\[ \left \{\frac {1}{\left (\frac {d}{d x}y \left (x \right )\right )^{7}}, \left (\frac {d}{d x}y \left (x \right )\right )^{2}, \left (\frac {d}{d x}y \left (x \right )\right )^{3}, \sqrt {\frac {d}{d x}y \left (x \right )}, \frac {d}{d x}y \left (x \right ) \right \} \]
Use indets with type identical(diff(y(x),x))^anything is used. But must use the flat
option to work correctly.
restart; expr:= y(x)*diff(y(x),x)^(1/3)+sin(x)*diff(y(x),x)^3 + z*diff(y(x),x)*ln(y(x)*diff(y(x),x)^2) +diff(y(x),x)^(1/2) + x/diff(y(x),x)^7; t1:=identical(diff(y(x),x))^anything; t2:=identical(diff(y(x),x)); indets[flat](expr, 'Or'( t1, t2 ));
gives
{diff(y(x), x)^(1/3), 1/diff(y(x), x)^7, diff(y(x), x)^2, diff(y(x), x)^3, diff(y(x), x), sqrt(diff(y(x), x))}
Without using flat it will given wrong result. For example
restart; expr:=diff(y(x),x)^2; t1:=identical(diff(y(x),x))^anything; t2:=identical(diff(y(x),x)); indets(expr, 'Or'( t1, t2 ));
Gives
{diff(y(x), x)^2, diff(y(x), x)}
You see, it has extra diff(y(x), x) showing up. Adding flat it gives
restart; expr:=diff(y(x),x)^2; t1:=identical(diff(y(x),x))^anything; t2:=identical(diff(y(x),x)); indets[flat](expr, 'Or'( t1, t2 ));
Now it gives
{diff(y(x), x)^2}
Which is the correct result.
To go from \(\ln (A B)\) to \(\ln A + \ln B\) need to use simplify with ln option but add assumptions that one of the
terms is positive. Else it will not do it
restart; simplify(ln(A*B),ln); # no change simplify(ln(A*B),ln) assuming A>0; # ln(A) + ln(B) simplify(ln(A*B),ln) assuming B>0; # ln(A) + ln(B)
To go from \(\ln ( \frac {A}{B})\) to \(\ln A - \ln B\) need to use simplify with ln option but add assumptions that
\(B>0\).
restart; simplify(ln(A/B),ln); # no change simplify(ln(A/B),ln) assuming A>0; # ln(A) + ln(1/B) simplify(ln(A/B),ln) assuming B>0; # do this: -ln(B) + ln(A)
To go from \(\ln A + \ln B\) to \(\ln (A B)\) need to use combine with assumptions that either \(A\) or \(B\) is positive, else it will
not do it.
restart; combine( ln(A) + ln(B),ln); # no change combine( ln(A) + ln(B),ln) assuming A>0; # ln(A*B) combine( ln(A) + ln(B),ln) assuming B>0 # ln(A*B)
To go from \(\ln A - \ln B\) to \(\ln \frac {A}{B} \) need to use combine with assumptions that either \(B\) is positive.
restart; combine( ln(A) - ln(B),ln); #no change combine( ln(A) - ln(B),ln) assuming A>0; # -ln(B/A) combine( ln(A) - ln(B),ln) assuming B>0 # use this ln(A/B)
Given list such as [1,2,3,4,5,x,y,8,9,Pi] find position of elements that are not numeric.
In this case the answer should be [6,7,10]
restart; lis:=[1,2,3,4,5,x,y,8,9,Pi]; lis2:=select(x->not(type(x,numeric)),lis); map(x->ListTools:-Search(x,lis),lis2) [6,7,10]
I could not find a way to do it using one command like with Mathematica. The first
command above uses select to first find non numeric entries. The second command
ListTools:-Search then find the index/position.
Maple’s ListTools:-Search should really have a version that allows one to select the
element directly. Something like this
lis:=[1,2,3,4,5,x,y,8,9,Pi]; ListTools:-Search(x->not(type(x,'numeric')),lis,all)
I needed to do this as I was translating Maple code to Sagemath. Where sagemath supports \(\arctan \) with only one argument.
Given an expression such as \[ x +\arctan \left (y,x \right )+\sin \left (x \right ) \] Convert it to \[ x +\arctan \left (\frac {y}{x}\right )+\sin \left (x \right ) \]
evalindets(expr,'specfunc(arctan)',f->`if`(nops(f)=2,arctan(op(1,f)/op(2,f)),f))
Given an expression such as \(a+\sin (x) + \pi + y+ f(r) \) how to find all symbols in it, which will be \(a,x,y,r\)?
One way
expr:=a+sin(x) + Pi + y+ f(r); indets(expr,And(symbol,Not(constant)));
Another is
expr:=a+sin(x) + Pi + y+ f(r); indets(expr,assignable(name));
Both give {a, r, x, y}
This question came up in https://mathematica.stackexchange.com/questions/274535/replacing-only-variables-in-specific-locations-with-replace-all
where the user wanted to I would like to replace R0 with r, but only in the arguments of functions.
The input they had is line = R0*f[R0,x] + R0^2*42*D[g[R0,x],x]
In Maple, this can be done as follows, which I think is a easier thanks to Maple’s strong type system.
line := R0*f(R0,x) + R0^2*42*diff(g(R0,x),x); evalindets(line, 'patfunc(identical(R0),anything)', Z-> subsop(1 = r, Z ));
Which gives
\[ \mathit {R0} f \left (r , x\right )+42 \mathit {R0}^{2} \left (\frac {\partial }{\partial x}g \left (r , x\right )\right ) \]
Given an expression which contains some different functions each with different number of arguments. Suppose we want to change only the last argument of each function if the last argument is \(x\), and change it to say \(x^2\).
how to do that?
Hence given R0*f(R0,x)+42*R0^2*D[2](g)(R0,x)+h(x,y,z,r,x) we want to change it to
R0*f(R0,x^2)+42*R0^2*D[2](g)(R0,x^2)+h(x,y,z,r,x^2).
expr:= R0*f(R0,x) + R0^2*42*diff(g(R0,x),x)+ h(x,y,z,r,x); evalindets(expr, 'patfunc[reverse](anything,identical(x))', Z-> subsop(-1 = x^2, Z ));
Note the use of subsop(-1 = x^2, Z ) where \(-1\) means the last entry in the argument of the
function. This is basically same as last example, but uses patfunc[reverse] instead of just
patfunc
This is the same example as the above, but now we want to remove the last argument instead of changing it.
Hence, given an expression which contains some different functions each with different number of arguments. Suppose we want to remove only the last argument of each function if the last argument is \(x\).
how to do that?
Hence given R0*f(R0,x)+h(x,y,z,r,x) we want to change it to R0*f(R0)+h(x,y,z,r).
expr:=R0*f(R0,x)+h(x,y,z,r,x) evalindets(expr, 'patfunc[reverse](anything,identical(x))', Z-> subsop(-1 = NULL, Z ));
Be careful using the above on expression that has diff(f(y,x),x) as this will give 0,
because we basically removed the variable of the differentiation.
Given \(3 \,{\mathrm e}^{x}+\sin \left (a \,{\mathrm e}^{x}\right ) f\left ({\mathrm e}^{5 x}\right )\) how to find all terms with pattern anything*exp ? If we do this
expr:=3*exp(x)+sin(a*exp(x))*f(exp(5*x)); indets(expr,`&*`(anything,'specfunc(exp)')); # {a*exp(x), 3*exp(x)}
We see it did not find exp(5*x) this is because there is nothing multiplying the exp
function. To find this we add Or like this to count for both cases
expr:=3*exp(x)+sin(a*exp(x))*f(exp(5*x)); indets[flat](expr,Or( `&*`(anything,'specfunc(exp)'), 'specfunc(exp)' )) # {a*exp(x), 3*exp(x), exp(5*x)}
The flat option is needed, as without it this will be the result
expr:=3*exp(x)+sin(a*exp(x))*f(exp(5*x)); indets(expr,Or( `&*`(anything,'specfunc(exp)'), 'specfunc(exp)' )) # {a*exp(x), 3*exp(x), exp(x), exp(5*x)}
Given an inert Sum such as \[ r = \sum _{n=2}^{\infty } a_{n -1} x^{n -1} \] How to obtain the body of the sum, the index variable, the lower starting value and the upper limit?
r:=Sum(a[n - 1]*x^(n - 1), n = 2 .. infinity); op(0,r) #head Sum op(1,r) #body a[n - 1] x^(n - 1) op(2,r) # sum specs n = 2 .. infinity lhs(op(2,r)) #name of summation index n rhs(op(2,r)) # lower..upper limits 2 .. infinity op(1,rhs(op(2,r))) #lower limit 2 op(2,rhs(op(2,r))) #upper limit infinity
Given a 1D container, such as vector or list, called \(A\), how to find if this sequence is inside another sequence say \(B\) and the indecies in \(B\) where the sequence \(A\) is located?
Use ArrayTool in version 2023
A:=[1,3,5]; B:=[1,3,4,5]; status,indices_list := ArrayTools:-IsSubsequence( A, B ,'output' = ['check','indices'],'match'='exact'); if status then print("Sequence ",A," Was found in ",B," At indices ",indices_list); else print("Sequence ",A," Was not found in ",B); fi;
The above gives
"Sequence ", [1, 3, 5], " Was found in ", [1, 3, 4, 5], " At indices ", [1, 2, 4]
If \(A\) was not sequence inside \(B\), then status will be false otherwise.
Use membertype
restart; expr:=cos(x)+x+9+sin(x)*(x^2+4); membertype( polynom(anything,x), expr,'loc' ); true loc; 2
The above says there is a polynomial in \(x\) inside the expression at op(2,expr) notice that \(x^2+4\) is
not a polynomial since the expression will expand and \(\sin x\) will be multiplied by it causing
it not to become polynomial. So only \(9+x\) is the polynomial. The location is where
the memeber starts at. Notice that Maple sorts polynomial from lower to higher
powers.
restart; expr:=cos(x)+x+9+sin(x)*(x^2+4); membertype( integer, expr,'loc' ); true loc; 3
The above says that there is an integer (which is \(9\) in this example) inside the expression at
op(3,expr)
I had case where I wanted to substitute _F=(y-x^2)/x in an expression obtained which is
x*diff(F1, x) + 2*y*diff(F1, y)
Using eval does not work
restart; eval(x*diff(F1, x) + 2*y*diff(F1, y),F1=(y-x^2)/x); 0
Also using delayed does not work
restart; eval('x*diff(F1, x) + 2*y*diff(F1, y)',F1=(y-x^2)/x); value(%) 0
ALso using subs does not work
restart; subs(F1=(y-x^2)/x,x*diff(F1, x) + 2*y*diff(F1, y)) 0
However, delayed with subs finally worked
restart; subs(F1=(y-x^2)/x,'x*diff(F1, x) + 2*y*diff(F1, y)'); value(%) x*(-2 - (-x^2 + y)/x^2) + 2*y/x
So the rule is, if you want to replace a function inside diff, use subs and not eval, and make
sure to delay evaluation of the expression and then use value() to obtain the final
result.
These are examples how to match expression types
expr is product a*b which is matched using
type(a*b,`&*`(anything,anything)) true
The above can also be written using infix notation
type(a*b, anything &* anything) true
However, the above only match product of two terms. To match 3*a*b use
type(3*a*b,`*`) true
expr is division a/b which is matched using same as `*`
type(a/b,`&*`(anything,anything)) true
Match type of a*f(x), which is anything times a function that takes one argument
\(x\).
type(a*f(x),`&*`(anything,patfunc(identical(x),anything))) true
Using patfunc is better than using function
type(a*f(x,y),`&*`(anything,function(identical(x)))) true
Because patfunc matches on \(f(x,y,z,\dots )\) and not just function which takes only one
argument \(f(x)\). But if you know your function takes only one argument, then use
function
To match 3*y/x. This was tricky. Had to use
expr:=3*y/x; type(expr, `&*`(anything,identical(y),`^`(identical(x),-1))) or type(expr, `&*`(identical(y),`^`(identical(x),-1))) true
I used or to account for possible expr:=y/x, i.e. missing constant at front. Note
that
expr:=3*y/x; type(3*y/x, `&*`(anything,identical(y/x))) false
Does not match, since 3*y/x is actually 3*y times 1/x internally.
Match on f(b* y/x)
type(f(3*y/x),function(`&*`(anything,identical(y),`^`(identical(x),-1)))) true
Match on f(b* y/x) or f(y/x)
expr:=f(3*y/x); type(expr, function(`&*`(anything,identical(y),`^`(identical(x),-1)))) or type(expr, function(`&*`(identical(y),`^`(identical(x),-1)))) true
Again, had to use or to account for missing constant multipler.
Match on a*f(b* y/x)
expr:=a*f(3*y/x); type(expr, `&*`(anything,patfunc(`&*`(anything,identical(y),`^`(identical(x),-1))))) true
Match against 9+3 y/x
expr:=9+3*y/x; selector:=`&+`(anything,`&*`(anything,identical(y),`^`(identical(x),-1))): type(expr, selector); true select(type,[expr],selector); [9 + 3*y/x]
Notice in the above, when using select we need to put the expression inside a list, as
select looks at each operand, This way the whole expression is taken as one.
If we just used select(type,expr,selector); it would not have found
it.
To use patmatch the command becomes
patmatch(expr, a::anything+b::anything*y/x,'la'); la [a = 9, b = 3]
The nice thing about patmatch is that it allowed one to assign variable to parts of the expression automatically.
Match against 9+f(3 y/x) where now \(f\) is function. Using patmatch. I could not do this
in one command, as all my attempts failed:
expr:=3+4*x*f(3*y/x); body_of_function:=C::anything*y/x; patmatch(expr,A::anything+B::anything*F::function(C::anything*y/x),'la'); patmatch(expr,A::anything+B::anything*F::patfunc(C::anything*y/x),'la'); Error, (in PatternMatching:-AlgStruct:-Match) testing against an invalid type Error, (in type/patfunc) testing against an invalid type
So I had to do it in two steps. First match on the function as whole, then use that to
match on f(3*y/x) in second stage, like this
expr:=3+4*x*f(3*y/x); patmatch(expr,A::anything+B::anything*F::function(anything),'la'); la; [A = 3, B = 4*x, F = f(3*y/x)]
And now
assign(la); A:='A'; patmatch (op(1,F),A::anything*y/x,'la'); la [A = 3]
Overall, I find Mathematica’s pattern matching constructs much simpler and more intutive to use and easier to learn as there are many examples and tutorials. For example, the last example above in Mathematica could be done as follows
expr=9+(4*x)*f[3*y/x]; Cases[{expr},any0_.+ any1_.*any2_[any3_.*y/x]:>{any0,any1,any2,any3}] {{9, 4 x, f, 3}}
Maple’s help pages are not good at all and provide little or no examples to learn from compared to Mathematica’s excellent help pages. For any serious pattern matching tasks, I would use Mathematica. Maple has a better debugger and hence easier to debug the code because of this. So it is a tradeff between these two systems.
I was trying to match on anyfunction that has \(x\) inside its arguments. It turned out that
matching \(f(a x)\) vs. \(f(a+x)\) needed to have the identical(x) being placed first when it is a sum and last
when it is a product. Very strange. Just be aware of this
indets(f(a+x),function(`&+`(anything,identical(x)))); #failed {}
But
indets(f(a+x),function(`&+`(identical(x),anything))); {f(x + a)}
While with product, it is the other way around
indets(f(a*x),function(`&*`(anything,identical(x)))); {f(a*x)}
But now this fail
indets(f(a*x),function(`&*`(identical(x),anything))); #fail {}
I have not yet figure how to tell it that the order does not matter. Maple 2023.1
A better way than the above, if I want to find any function that takes in \(x\) or \(y\) as arguments is
to use patfunc like this
indets(f(a*x+b+y),patfunc(anything,`Or`({identical(x),identical(y)}))); {f(a*b + y)} indets(f(b+y),patfunc(anything,`Or`({identical(x),identical(y)}))); {f(b + y)} indets(exp(a*x),patfunc(anything,`Or`({identical(x),identical(y)}))); {exp(a*x)} indets(exp(a+x),patfunc(anything,`Or`({identical(x),identical(y)}))); {exp(x + a)}
restart; with(MmaTranslator); #load the package FromMma(`Integrate[Cos[x],x]`);
Or
restart; with(MmaTranslator); #load the package convert(`Integrate[Cos[x],x]`, FromMma);
f:=proc() eq:=x*diff(y(x),x)+y(x)=exp(2*x); dsolve(eq,y(x)); end proc;
Then used the command stopat(f); then called the procedure f(); and now the debugger
comes up. Did step command and now it steps inside dsolve
Some examples
stopat(`ODEtools/symtest`); stopat(`ODEtools/test`); stopat(`ODEtools/normal/expanded`); stopat(`ODEtools/odepde`); stopat(`ODEtools/odeadv`); #for DEtools:-odeadvisor stopat(`odsolve/dAlembert`); stopat(`odsolve/dAlembert/integrate`); stopat(`odsolve/answer`); stopat(`odsolve/homogeneous`); #for all A,C,D,G types stopat(`odsolve/homogeneous_C/integrate`); stopat(`odsolve/exact`); #for solving exact ODE stopat(`odsolve/exact/integrate`); stopat(`odsolve/exact/integrate`(f,y,x,M,N)); #where f here is RHS of y'=RHS; DEtools:-symtest([-3,y],ode,y(x)); `ODEtools/normal/expanded`
For exact ode, can also do
ode:=....# write your ode here Student:-ODEs:-Solve:-Exact(ode,y(x),output=steps);
For integration use
infolevel[`evalf/int`]:=5;infolevel[int]:=5;
Another option
restart; interface(verboseproc=3) #(try 2 also)
then print(procedure); or eval(procedure_name); for example
restart: interface(verboseproc=3): print(LinearAlgebra:-GramSchmidt); print(lcm);
Also can use showstat, in this case interface(verboseproc=3) is not needed. Also
showstat gives line numbers and I think it is easier to read. Some examples
showstat(`odsolve/2nd_order`) showstat(`evalf/hypergeom`); showstat(`evalf/exp/general`); showstat(`evalf/Psi`); showstat(`evalf/int`); showstat(`dsolve/SERIES`); #these 3 shows the main 3 functions by each solver showstat(`odeadv/dAlembert`); #used by advisor showstat(`odsolve/dAlembert`); # main API. showstat(`odsolve/dAlembert/integrate`); #used to integrate the ode showstat(`ODEtools/odeadv`); showstat(DEtools:-odeadvisor); showstat(`dsolve/series/froben/inhom`) showstat(`dsolve/series/froben`)
To stop at anyone of these functions in debugger do
stopat(`dsolve/series/froben/inhom`) #code here, say dsolve command.
The above will stop in the debugger in the above function.
There is also a function by Joe Riel here here is the post by Joe Riel:
"A disadvantage of showstat, particularly if you want to cut and paste the output, is that it includes line numbers. Here is a simple procedure I threw together to remove the line numbers."
PrintProc := proc(p::name,lines::{posint,posint..posint}) local width; option `Copyright (C) 2004 by Joseph S. Riel. All rights reserved.`; description "Print like showstat, but without line numbers"; width := interface('screenwidth'=200); try printf("%s", StringTools:-RegSubs( "\n ...." = "\n" ,debugopts('procdump'= `if`(nargs=1,p,[args])))) catch "procedure name expected": error "%1 is not a procedure name",p finally interface('screenwidth'=width) end try; NULL end:
To print source code to file using the above, do the following
currentdir("C:\\data"); interface('prettyprint'=1): interface('verboseproc'=3): writeto("listing.txt") PrintProc('singular'); writeto('terminal'):
Now the output will show up in the file "listing.txt" and also no line wrapping. The above I found is the best solution so far to do this.
trace(foo); untrace(foo);
also see debug(foo);
Also
infolevel[all]:=5: printlevel:=10:
See http://www.mapleprimes.com/questions/35951-How-To-Debugtrace-Things-In-Maple
Also look at kernelopts(opaquemodules=true)
Here is a useful post by Carl Love from Maple prime forum that summarizes all of these
Here are four things that you can do to get more information. I have listed them in order by how structured the information is, with the most structured first.
Set
infolevel[all]:= 5;
That will cause programs to print out additional information of the programmers’ choosing. You can use higher or lower numbers for more or less information. Most programs don’t use levels higher than 5.
Print the code of procedures with showstat:
showstat(int); showstat(sin); showstat(cos);
Trace the execution of particular procedures with trace:
trace(int); trace(sin);
Trace the execution of everything with printlevel:
printlevel:= 10000:
You can use higher or lower numbers for more or less information.
Some examples
interface(verboseproc=3); print(DEtools) print(`ODEtools/symgen`); print(`symgen/methods`); print(`symgen/do`);
To stop the debugger at symgen do
stopat(`ODEtools/symgen`);
To get infolevel on symgen do
infolevel[`symgen`]:=5;
Or to see line numbers
interface(verboseproc=3); showstat(dsolve)
Or can use the Browse(); command
with(LibraryTools); Browse();
Another option I found is
s:=debugopts(procdump=`showstat`);
Then the above produces listing that can be copied as string with line wrapping ok.
One way
L:=[]: for i from 1 to 3 do : L:=[op(L),i]; end do;
But a better way is to use seq if one knows the length
L:=[seq(i,i=1..3)]; L := [1, 2, 3]
Since list is unmutable, a more efficient method, for long lists, is to use Array, and then convert the result back to list at the end since Array can grow dynamically without preallocation each time something is inserted as follows
L:=Array(): for i from 1 to 3 do : L(i):=i; end do; for i from 1 to numelems(L) do : print(L[i]); end do; L := convert(L,list)
Which wil print
L := [1] L := [1, 2] L := [1, 2, 3] 1 2 3 L := [1, 2, 3]
Notice that to add to an Array, () is used. But to access an entry in an array [] is
used.
And finally, using Array also, it can be done without using any indexing as follows
L:=Array(1..0): for i from 1 to 3 do : L ,= i; end do; L := convert(L,list)
For the above to work, the array must be declared using Array(1..0). The new syntax
A ,= i will append to the array, and there is no need to write A(i) := i
By Carol Devore on the net:
Use infolevel. For example, to show what logic dsolve uses, do this: First try > infolevel[all]:= 5; That will probably give more information than you want, but if not, then try > printlevel:= 1000; If you want information about a specific procedure, you can use debug. For example, restart; debug(`int/int`); int(p, x= 0..1); To find out what procedures are being called without getting too much extra information, use excallgraph.
Trying on dsolve
infolevel[dsolve]:= 3; dsolve({eq1},y(x)); Methods for second order ODEs: Trying to isolate the derivative d^2y/dx^2... Successful isolation of d^2y/dx^2 --- Trying classification methods --- trying a quadrature trying high order exact linear fully integrable trying differential order: 2; linear nonhomogeneous with symmetry [0,1] trying a double symmetry of the form [xi=0, eta=F(x)] <- double symmetry of the form [xi=0, eta=F(x)] successful
Here, I am looking at fouries series expansion of \(f(x)=0\) between \(–\pi \) and 0, and \(f(x)=1\) between 0 and \(\pi \).
The Fouries series expansion is worked out to be as below. This shows that the series approximate the above \(f(x)\) as more terms are added
restart; f:=(x)-> 1/2 + (1/Pi)*(sin(x)+sin(3*x)/3+sin(5*x)/5+sin(7*x)/7); plot(f(x),x=-10..10);
From DOS, point to where your cmaple is
>"C:\Program Files\Maple 7\BIN.WNT\"cmaple
To make it execute maple commands use the < foo.txt to pipe maple commands in the file
to it.
A:= Matrix( [ [1, 2, 3] , [3, 6, 7] , [5, 6, 9] , [7, 7, 7] ]); whattype(A); Matrix size:=LinearAlgebra:-Dimension(A); size := 4, 3 row:=size[1]; row := 4 col:=size[2]; col := 3
You can extract any part of the matrix like this:
B:=A[1..3,2..2];
\[ \left [ \begin {array}{c} 2\\ 6\\ 6\end {array} \right ] \]
By Carl Devore http://mathforum.org/kb/message.jspa?messageID=1570678
Maple list and sequence structures are more flexible than Matrices, which are highly structured. A Maple list of lists (called a listlist in Maplese) is akin to a matrix in some other languages. Many matrix operations can be performed directly on the listlist form, but to do serious linear algebra, you should convert to a Matrix. Of course, it is trivial to convert a listlist to Matrix: LL:= [[1,2], [3,4]]; M:= Matrix(LL); So here is another solution in line with your original wishes. This is "index free", but the table-based solution I gave earlier should be faster. (It is usually considered bad form to repeatedly append to a list or sequence.) L:= [][]; # Create a NULL sequence do line:= readline(file); if line::string then if line contains valid data then Z:= a list of that data; L:= L, Z fi else break fi od A:= Matrix([L]); # Note []: seq -> list.
To move move a column into a matrix: Here, I want to copy 2nd column to the 3rd column:
A; \[ \left [ \begin {array}{ccc} 1&2&3\\ 3&6&7\\ 5&6&9\\ 7&7&7 \end {array} \right ] \]
B:=A[1..row,2]; \[ \left [ \begin {array}{c} 2\\ 6\\ 6\\ 7 \end {array} \right ] \]
A[1..row,3]:=B: A;
\[ \left [ \begin {array}{ccc} 1&2&2\\ 3&6&6\\ 5&6&6\\ 7&7&7 \end {array} \right ] \]
Maple can return multiple values. Make sure to use the comma "," in the body of the procedure to separate each return value. Example:
size_matrix:=proc(x) 3*x, 4*x; end proc; row,col :=size_matrix(5);
When passing a variable to maple procesure, the variable VALUE is passed to the procedure (This is different from say Fortran where the default is pass by reference). But this is the same as with Mathematica.
For example, if a variable X had value 10, then you call a procedure FOO passing
it X, then inside FOO, X will be the number 10, not the argument variable X.
So, this means one can not have X on the left hand side inside FOO. Like this
x:=1
The only way to assign new value to the input and return new value, is to use a local variable, like this:
one:= proc(x) local y; print(x); y:=x+ 1; print(x); y; end proc; z:='z'; z:=5; f:=one(z); f := 6
Use `type/name` to define new type name.
`type/char`:= x-> x::string and length(x)=1; P:= proc(c::char) print(c) end proc: P("x"); "x" P("xy"); Error, invalid input: P expects its 1st argument, c, to be of type char, but received xy > `type/byte`:= x-> x::integer and (x>= 0 and x<256); #will define a byte (unsigned integer)
Code from net by Carl Devore:
MMax:= proc(M::{Matrix,matrix}) local C,r,c,mx,L,p; C:= op(`if`(M::Matrix, [1,2], [2,2,2]), eval(M)); L:= map(op, convert(M, listlist)); mx:= max(L[]); member(mx,L,'p'); r:= iquo(p, C, 'c'); mx, `if`(c=0, [r,C], [r+1,c]) end;
Code below from C W
A:=matrix(12,12,rand(100)); Ao:=array((proc(E) local i; [seq(i=(rhs=lhs)(E[i]),i=1..nops(E))]end) (sort(op(3,eval(A)),proc(E1,E2) if rhs(E1)>rhs(E2) then true else false fi end))); Ao[1];
First create the module:
restart; nma:= module() option package; export getMaxMatrix; getMaxMatrix := proc (M::{matrix, Matrix}) local C, r, c, mx, L, p; C := op(`if`(M::Matrix,[1, 2],[2,2,2]),eval(M)); L := map(op,convert(M,listlist)); mx := max(L[]); member(mx,L,'p'); r := iquo(p,C,'c'); mx, `if`(c = 0,[r, C],[r+1, c]) end proc; end module; A:= Matrix( [ [1, 2, 3] , [3, 6, 7] , [5, 6, 9] , [7, 7, 7] ]); nma[getMaxMatrix](A);|
Gives 9, [3, 3]. Now save the module.
savelibname := "C:/MAPLE_PACKAES"; march('create', savelibname, 20);
now save the library to disk. savelib(nma);
Now we can test everything by reinitialize everything and reload the library.
>restart #Add my library to LIBNAME >libname:="C:/MAPLE_PACKAGES",libname; > A:=matrix( [ [1,2,3],[4,6,9] ]); >with(nma); >nma[getMaxMatrix](A);
Now to print a proc() in the package, do
>interface(verboseproc=3); > print(nma[getMaxMatrix]);
Now you can list what packages exist in the archive:
march('list',savelibname); march('extract',savelibname,":-1.m","C:MAPLE_PACKAGES/t.m")
Some notes. need to clean later
> module1lib:=`module1\\lib`; > system("md "||module1lib); > march('create',module1lib,100); > makehelp(module1,`module1/module1.mws`,module1lib): > makehelp(`module1/export1`,`module1/export1.mws`,module1lib): > savelibname:=module1lib: ### doesn't affect current libname > savelib(module1); ### no error message > restart; > module1lib:="module1\\lib": > libname:=module1lib,libname; ### now Maple will find module1 > with(module1); > ?module1
Also there is a long thread here on Maple prime on making personal packages in Maple How-To-Create-A-Personal-Package
From: Robert Israel (israel@math.ubc.ca) Subject: Re: Getting non-integral results in hex Newsgroups: comp.soft-sys.math.maple Date: 2003-06-13 00:07:37 PST I assume you mean floating-point numbers. Note that Maple floats (as opposed to "hardware floats") are in fact stored in base 10. To convert a float to hex with n digits after the ".", you can use this: > `convert/hexfloat`:= proc(x::numeric, n::nonnegint) local A,B,ax,R; if nargs = 1 then return procname(x,round(Digits*log[16](10))) fi; if x = 0 then return cat(`0.`,`0`$n) fi; ax:= abs(x); A:= floor(ax); B:= round(frac(ax)*16^n); if B = 16^n then A:= A+1; B:= 0 fi; R:= cat(convert(A,hex),`.`); if x < 0 then R:= cat(`-`,R) fi; cat(R,substring(convert(16^n+B,hex),2..-1)); end; And then, e.g.: > convert(1234.5678, hexfloat, 4); 4D2.915B
mtaylor(sin(x),[x],10);
\[ x-1/6\,{x}^{3}+{\frac {{x}^{5}}{120}}-{\frac {{x}^{7}}{5040}}+{\frac { {x}^{9}}{362880}} \]
restart; a:=Matrix([ [2,3,4],[4,5,6] ]); nRow,nCol :=LinearAlgebra[Dimension](a); for i from 1 to nRow do for j from 1 to nCol do printf("a(%d,%d)=%d\n",i,j,a[i,j]); end do; end do; a(1,1)=2 a(1,2)=3 a(1,3)=4 a(2,1)=4 a(2,2)=5 a(2,3)=6
restart; a:=Matrix([ [2,4],[5,7] ]); LinearAlgebra:-Determinant(a); -6
H := LinearAlgebra:-HilbertMatrix(5);
\[ \left [ \begin {array}{ccccc} 1&1/2&1/3&1/4&1/5\\ \noalign {\medskip }1/ 2&1/3&1/4&1/5&1/6\\ \noalign {\medskip }1/3&1/4&1/5&1/6&1/7 \\ \noalign {\medskip }1/4&1/5&1/6&1/7&1/8\\ \noalign {\medskip }1/5&1/6&1 /7&1/8&1/9\end {array} \right ] \]
Matlab is much easier here. In maple, need to covert the matrix to a list of list of points first.
restart; H := LinearAlgebra:-HilbertMatrix(5): nRow,nCol :=LinearAlgebra[Dimension](H): L:=[seq([seq( [i,j,H[i,j]], i=1..nRow) ], j=1..nCol)]: plots:-surfdata(L);
An error in maple raises an exception. So, use try catch to trap it as follows:
try v,pos:=MMax(4); catch: printf("an error is cought\n"); end try;
From the net, by Carl Devor:
`print/commas`:= proc(N::integer) local n,s,i,b; n:= ListTools:-Reverse(convert(abs(N), base, 1000)); if N<0 then n:= subsop(1= -n[1], n) fi; nprintf("%s", sprintf(cat("%d", ",%03d" $ nops(n)-1), n[])) end proc: commas(456554); 456,554
To convert a string to array of chars use array(StringTools:-Explode(S))
s:="Nasser M. Abbasi": r:=array(StringTools:-Explode(s)); r:=["N" "a" "s" .......]
Now can use the string as normal array
r[4]; "s"
Units[GetDimensions](base); amount_of_information, amount_of_substance, currency, electric_current, length, logarithmic_gain, luminous_intensity, mass, thermodynamic_temperature, time
Use the Sum command.
restart; expr:= (-1)^i/(2*i+1)^2; Sum(expr,i=0..infinity); evalf(%,50); 0.91596559417721901505460351493238411077414937428167
Notice, if I used the sum command instead of the Sum command I get this result:
sum(expr,i=0..infinity); Catalan
This shows how to do a simple package and use it without building a library. Just using a plain text file.
Create this nma_pkg1.txt file:
nma_pkg1 := module() export f1; option package; f1:= proc() print("in pakcage nma_pkg1"); end proc; end module;
now save it, and from maple do
>read("c:\\nma_pkg1.txt");
now execute f1() as this:
>nma_pkg1[f1](); "in pakcage nma_pkg1"
now put it in a library (so that we can use with, instead of read)
> savelibname:=("c:/maple"); > march('create', savelibname, 20); > savelib(nma_pkg1); >restart; > libname := "c:/maple",libname; > with(nma_pkg1); > f1(); "in pakcage nma_pkg1"
now make changes to the nma_pkg1.txt file and updated again as above.
?index,package
restart; f:=3*x^2 + y* cos(x*y); the_grad :=linalg[grad](f,[x,y]); plots[fieldplot](the_grad,x=-2..2,y=-2..2);
or
or can do it in just one command: plots[gradplot](f,x=-2..2,y=-2..2);
Suppose you want the 100 digits of Pi put in a list. This is one way to do it:
restart; L:=evalf(Pi,100); S:=convert(L,string); the_list:=[seq(parse(S[i]),i=3..length(S))]; the_list := [1, 4, 1, 5, 9, 2, 6, 5, 3, ..
This below now tells how many times each digits occurs.
>stats[transform,tally](the_list); [Weight(0, 8), Weight(1, 8), Weight(2, 12), Weight(3, 11), Weight(4, 10), Weight(5, 8), Weight(6, 9), Weight(7, 7), Weight(8, 13), Weight(9, 13)]
Written sometime in 2005? I should really record the time when I write something.
I just run these now, Auust 2014, and now Maple 18 as very fast. So this all below is no longer valid. I will leave it here for now for reference until I update it all later
I have written a few lines of code, which counts how many times each digit occurs after the decimal points of \(\pi \)
Written this in maple first. Then did similar thin in mma 5.0. Both are run on the same PC. No other applications are running at the time when I run the code.
The basic idea of the algorithm is to use evalf(Pi,digits) in maple to find \(\pi \) for any
number of decimal digits, and to use N[Pi,digits] in mma for doing the same. (Where the
variable digits above is the number of digits)
Then in maple convert the above \(\pi \) to a string, and generate a sequence of the characters
to right of decimal point, then use stats[transform,tally] to do the actual
counting.
In mma, I use RealDigits[] to get a list of the digits, and then use Count[] to do the
counting.
This is result of some of the runs to find Pi to some digits, and the total time (to find Pi and do the counting)
All times are in cpu seconds, machine is P4, 2.8 Ghz, 500 MB of RAM, single CPU, hyperthreading enabled, running XP home edition. Maple 9.03 student version, and mma 5.0 student version.
Below is the result, and below that I show the maple code and the mma code.
Because of this, before each run in mma, I exited the application and started it fresh. In maple, it does not matter for the above reason.
100,000 digits: Find_Pi Total Maple 9.0 55 84 Mma 5.0 0.9 1.54
Mma is 60 times faster in finding pi and about 56 times faster overall
300,000 digits: Find_Pi Total Maple 9.0 309 781 Mma 5.0 3.7 6
Mma is 300 times faster in finding Pi, and 130 times faster overall.
3,000,000 digits Find_Pi Total Maple 9.0 Mma 5.0 85 118 Maple time in hours ! Still running.
Maple code
> restart; startingTime :=time(); L:=evalf(Pi,100000): timeToFindPiInSecs:=time()-startingTime; S:=convert(L,string): the_list:=[seq(parse(S[i]),i=3..length(S))]: stats[transform,tally](the_list); endingTime :=time(): cpuTimeInSecs := endingTime - startingTime;
mma code
Clear[] startingTime=TimeUsed[] t1=N[Pi,100000]; timeToFindPiInSecs=TimeUsed[]-startingTime {c,d}=RealDigits[t1]; theList=c[[Range[2,Length[c]]]]; f[digit_]:=Count[theList,digit]; r=Range[0,9]; Map[f,r] cpuTimeInSecs=TimeUsed[]-startingTime
update 12/25/03 Changed maple code on how to do the counting : To use
StringTools[CharacterFrequencies](S)
Now the counting in maple is much faster. It is always hard to know which is the best function to use.
restart; startingTime :=time(); L:=evalf(Pi,300000): timeToFindPiInSecs:=time()-startingTime; S:=convert(L,string): StringTools[CharacterFrequencies](S); endingTime :=time(): cpuTimeInSecs := endingTime - startingTime;
From: Ken Lin (maplemath@tp.edu.tw) Subject: Re: how to find which package a function belongs to? Newsgroups: comp.soft-sys.math.maple Date: 2003-12-04 03:49:26 PST When Maple first loaded, There are only two kinds of "internal" commands which can be called directly. One is the "kernal" commands coded in C, and the other includes many "internal" prodecures programmed by the kernal commands which lies in the "Main Library", There are also many other "external" procedures which were categorized into so called "packages", plots[display](...) for example, plots[] is a package(Library), and display() is the prodecure inside plots[]. All the packages can be loaded by with() command, like > with(plots); Because Different Packages include user library might have the same procedure name, Maple doesn't realize the "procedure_name" you type in, it took it for a "symbol". If you really want to know which packages provided by Maple the external procedure lies in, just mark the procedure_name and press F1 key, the Maple Help Browser will show you the packages you might be interested. By the way, plot3d() is a "internal" procedure lies in the Main Library. You can confirm that by: > op(0, eval(plot3d)); procedure or in Maple 9 > type( plot3d, 'std' ); #Is it internal? true > type( plot3d, 'stdlib' ); #Does is lie in "Standard(Main) Library"? true If you are interested the codes inside plot3d()... > interface(verboseproc=2): #Turn on verboseproc > print(plot3d); #eval() also works > interface(verboseproc=1): #Turn off verboseproc I hope this will give you some help. Have fun with Maple. Ken Lin
restart; f:= t->sin(omega*t) ; L:=convert(inttrans[laplace](f(t),t,s),int);
\[ {\frac {\omega }{{\omega }^{2}+{s}^{2}}} \]
To find the inverse, do:
inttrans[invlaplace](L,s,t);
\[ \sin \left ( \omega \,t \right ) \]
For unit step, use
_EnvUseHeavisideAsUnitStep:=true; f:=Heaviside(t-a); INV:=inttrans:-laplace(f,t,s) assuming a>0; #make sure to use a>0
\[ \frac {{\mathrm e}^{-s a}}{s} \]
Another example
_EnvUseHeavisideAsUnitStep:=true; f:=Heaviside(t)-Heaviside(t-a); INV:=inttrans:-laplace(f,t,s) assuming a>0; #make sure to use a>0
\[ \frac {1-{\mathrm e}^{-s a}}{s} \]
Any difference between using `diffalg/Rosenfeld_Groebner`(args) or diffalg[Rosenfeld_Groebner](args)
restart; f:= (x,y)->x^3-3*x*y^2; plot3d(f,-1..1,-1..1,numpoints=2500,style=patchcontour);
Use map
map(`^`,{1,2,3},3); {1, 8, 27}
incr:=.25; start:=0; last:=3; seq(start+i*incr,i=1..(last/incr));
read ?MVshortcut, ?MVassignment, and ?Mvextract and Transpose(R) can be shortened
to R^%T
Written feb 20, 2004
This is problem 7.4 chapter 4, in the Mary Boas book. Given \begin {align*} x s^2+y t^2 &= 1\\ x^2 s+y^2 t &= xy-4 \end {align*}
Find \(\frac {dx}{dt}, \frac {dx}{ds}, \frac {dy}{dt}, \frac {dy}{ds}\) at \(x=1,y=-3,s=2,t=-1\)
This is how I did it in maple:
restart; alias(x=x(s,t)); alias(y=y(s,t)); alias(Xt= diff(x(s,t), t)); alias(Xs= diff(x(s,t), s)); alias(Yt= diff(y(s,t), t)); alias(Ys= diff(y(s,t), s)); eq1:= x*s^2+y*t^2=1; eq2:= x^2*s+y^2*t=x*y-4; r1:=diff(eq1,t); r2:=diff(eq1,s); r3:=diff(eq2,t); r4:=diff(eq2,s); sol:=solve({r1,r2,r3,r4},{Xt,Xs,Yt,Ys});
\begin {align*} {\frac {\partial }{\partial s}}x \left ( s,t \right ) &= -{\frac {x \left ( s,t \right ) \left ( x \left ( s,t \right ) {t}^{2}-4\,y \left ( s,t \right ) st+2\,x \left ( s,t \right ) s \right ) }{2\,x \left ( s,t \right ) s{t}^{2}-2\,y \left ( s,t \right ) t{s}^{2}+x \left ( s,t \right ) {s}^{2}-y \left ( s,t \right ) {t}^{2}}}\\ {\frac {\partial }{\partial t}}x \left ( s,t \right ) &=-{\frac {y \left ( s,t \right ) t \left ( -3\,y \left ( s,t \right ) t+2\,x \left ( s,t \right ) \right ) }{2\,x \left ( s,t \right ) s{t}^{2}-2\,y \left ( s,t \right ) t{ s}^{2}+x \left ( s,t \right ) {s}^{2}-y \left ( s,t \right ) {t}^{2}}}\\ {\frac {\partial }{\partial s}}y \left ( s,t \right ) &=-{\frac {x \left ( s,t \right ) \left ( 3\,x \left ( s,t \right ) s-2\,y \left ( s,t \right ) \right ) s}{2\,x \left ( s,t \right ) s{t}^{2}-2\,y \left ( s,t \right ) t {s}^{2}+x \left ( s,t \right ) {s}^{2}-y \left ( s,t \right ) {t}^{2}}}\\ {\frac {\partial }{\partial t}}y \left ( s,t \right ) &=-{\frac {y \left ( s,t \right ) \left ( 4\,x \left ( s,t \right ) st-y \left ( s,t \right ) {s }^{2}-2\,y \left ( s,t \right ) t \right ) }{2\,x \left ( s,t \right ) s{t} ^{2}-2\,y \left ( s,t \right ) t{s}^{2}+x \left ( s,t \right ) {s}^{2}-y \left ( s,t \right ) {t}^{2}}} \end {align*}
points:= {x=1,y=-3,s=2,t=-1}; subs(points,sol);
This is problem 7.15 chapter 4 in Boas:
Given \(x^2 u-y^2 v=1\) and \(x+y=uv\) Find \(\frac {dx}{du},v\) and \(\frac {dx}{du},y\)
This is the maple code to solve this:
restart; eq1:=x^2*u-y^2*v=1; eq2:=x+y=u*v; r1:=D(eq1); r2:=D(eq2); r1_:=subs(D(v)=0,r1); r2_:=subs(D(v)=0,r2); sol:=solve({r1_,r2_},{D(x),D(u)}); print("dx/du,v="); rhs(sol[1])/rhs(sol[2]); r1_:=subs(D(y)=0,r1); r2_:=subs(D(y)=0,r2); sol:=solve({r1_,r2_},{D(x),D(u)}); print("dx/du,y="); rhs(sol[1])/rhs(sol[2]);
by http://www.math.fsu.edu/~bellenot
restart; t2 := proc(i, x, y) if i < 2 then [[x, y], [x, y - 1]], [[x, y], [x + 2^i, y - 1]] else [[x, y], [x, y - 1]], [[x, y], [x + 2^i, y - 1]], t2(i - 1, x, y - 1), t2(i - 1, x + 2^i, y - 1) end if end proc; PLOT(CURVES(t2(6,0,0)));
restart; z:= Int( sin(t)/t, t=sin(x)..cos(x)); diff(z,x);
\[ -{\frac {\sin \left ( x \right ) \sin \left ( \cos \left ( x \right ) \right ) }{\cos \left ( x \right ) }}-{\frac {\cos \left ( x \right ) \sin \left ( \sin \left ( x \right ) \right ) }{\sin \left ( x \right ) }} \]
restart; c:='c': C:='C': n:='n': P:='P': C := n -> ((n+2)/(3*n+1))^n: ### WARNING: calls to `C` for generating C code should be replaced by codegen[C] `The general term is `, c[n]= C(n); ` `; `The n-th root is:`; ### WARNING: calls to `C` for generating C code should be replaced by codegen[C] P := C(n)^(1/n): abs(c[n])^(1/n) = P; P := simplify(P, assume=positive): abs(c[n])^(1/n) = P;
restart; f:= 1/( (1-2*z)*(5*z-4) ); residue(f,z=4/5);
\[ \frac {-1}{3} \]
_EnvAllSolutions:=true; solve(sin(x)=0);
Pi _Z1~
Subject: Associated Legendre Author: Mehran Basti <Basti@worldnet.att.net> Organization: AT&T Worldnet Date: Mon, 25 Nov 2002 02:48:15 GMT Dear newsgroup: I had mentioned that my methods will solve classical equations without the use of infinite series. The following is a Maple code of my old files. Those days I had Maple2 but the general idea is the same in the process and you see that we can also solve the integrals involved. It does not make sense how are the theory behind it but eventually it will come into light. Just read the procedures and you can see the solution of associated legendre AL at the end. > s1:=-diff(p(t),t)+p(t)^2; > > s2:=exp(2*int(p(t),t))*T(t); > s3:=s1+s2; > s4:=diff(T(t),t)/T(t); > s5:=-(1/2)*(diff(s4,t))+(1/4)*s4^2; > s6:=s5+s2; > p(t):=-1/t+(1)/(2-t); > s1:=simplify(s1); > s1:=collect(%,t); > s2:=simplify(s2); > s1+s2=(2*t^2-4*t+m^2-1)/(t*(-2+t))^2; > solve(%,T(t)); > T(t):=simplify(%); > s2:=simplify(s2); > s2+s1; > s3:=simplify(%); > > s6:=simplify(s6); > t*(-2+t); > simplify(%); > z:=(r3*t^3+r2*t^2+r1*t+r0)/(%); > > simplify(diff(z,t)+z^2-s6); > s7:=collect(numer(%),t); > > coeff(%,t,0); > solve(%,r0); > r0:=op(1,{%}); > coeff(s7,t,1); > solve(%,r1); > r1:=simplify(%); > coeff(s7,t,2); > solve(%,r2); > r2:=simplify(%); > coeff(s7,t,3); > solve(%,r3); > r3:=simplify(%); > simplify(s7); > s3:=simplify(s3); > s4:=simplify(s4); > s6:=simplify(s6); > T(t):=simplify(T(t)); > z:=simplify(z); > 1/2*s4+2*p(t)+z; > s8:=simplify(%); > exp(int(%,t)); > expand(%); > g:=(%); > simplify(g,power); > g:=%; > Int(%,t); > Integralg:=(%); > int(g1(t),t); > x1:=-p(t)+g1(t)/(%); > diff(x1,t)+x1^2-s3; > simplify(%); > s10:=numer(%); > solve(%,int(g1(t),t)); > Ing:=(%); > simplify(subs(g1(t)=g,%)); > > Ing:=(%); > expand(%); > Ing:=simplify(%); > simplify(diff(%,t)-g); > expand(%); > simplify(%); > x:=-p(t)+g/Ing; > simplify(diff(x,t)+x^2-s3); > int(x,t); > exp(%); > expand(%); > s11:=simplify(%); > ALT:=t*(2-t)*diff(u(t),t$2)+2*(1-t)*diff(u(t),t)+(2-m^2/(1-(1-t)^2))*u(t); > -2*(1-t)/(2*t*(2-t)); > int(%,t); > exp(%); > s12:=simplify(%,power); > > u1:=s12*s11; > u1:=simplify(%,power); > simplify(subs(u(t)=u1,ALT)); > AL:=(1-nu^2)*diff(u(nu),nu$2)-2*nu*diff(u(nu),nu)+(2-m^2/(1-nu^2))*u(nu); > > u2:=subs(t=1-nu,u1); > simplify(subs(u(nu)=u2,AL)); > The advantage of these methods are that there are ample rooms for advances. Today my skills for solving classical equations such as Riccati is much advanced. Highly complicated and more general Riccati equations in its billions now possible. Sincerely Dr.M.Basti
To plot mapping of complex function in maple, use [plots]conformal The trick is to how
to specify the quadrant in the x-y plane. This example shows how.
Suppose we want to map the first quadrent. Then we specify the DIAGONAL points in the
range, from the lower left corner to the upper right corner, which then should be 0..1+I
Because 0 is the lower left corner, and \((1,i)\) is the upper right corner. Example:
restart; assume(y,real); assume(x,real); #f:= z->I+z*exp(I*Pi/4); f:= z->z^2; w:=f(x+I*y); u:=Re(w); v:=Im(w); plots:-conformal(f(z),z=0..1+I,grid=[16,16],numxy=[16,16],scaling=constrained);
This below uses the first TWO quadents, i.e. the upper half of the x-y plane
restart; assume(y,real); assume(x,real); #f:= z->I+z*exp(I*Pi/4); f:= z->z^2; w:=f(x+I*y); u:=Re(w); v:=Im(w); plots:-conformal(f(z),z=-1-I..1+I,grid=[16,16],numxy=[16,16],scaling=constrained);
This below puts the plots next to each others so to see them
restart; assume(y,real); assume(x,real); f:= z->I+z*exp(I*Pi/4); #f:= z->z^2; w:=f(x+I*y); u:=Re(w); v:=Im(w); A := array(1..2): A[1]:=plots:-conformal(z,z=0..1+I/2,grid=[16,16],numxy=[16,16],scaling=constrained): A[2]:=plots:-conformal(f(z),z=0..1+I/2,grid=[16,16],numxy=[16,16],scaling=constrained): plots:-display(A);
interface(showassumed=0) removes all tildas and interface(showassumed=1) adds the
tildas.
I wrote this to generate FS in Maple for some HW I was doing. I think this was for Math 121A at UC Berkeley in 2003
restart; f:=x->piecewise(-Pi<x and x<Pi/2,-1, Pi/2<x and x<1,0,1); assume(n,integer); nmaFourier2:=proc(f,freq,from_,to_,maxN) local n::integer,denomC,denomS,a,b; denomC:=( to_ - from_ ) / 2; denomS:=( to_ - from_ ) / 2; a:=proc(n) int(f(x)*cos(n*freq*x),x=from_..to_) /denomC; end proc; b:=proc(n) int(f(x)*sin(n*freq*x),x=from_..to_) / denomS; end proc; evalf(denomC); 1/2*a(0) + sum( a(n) * cos(n*freq*x) ,n=1..maxN) + sum( b(n) * sin(n*freq*x) ,n=1..maxN) end proc; r:=[seq(nmaFourier2(f,1,-Pi,Pi,nIter),nIter=1..10)]; plot(r,x=-Pi..Pi);
To animate do
g:=n->plot(nmaFourier2(f,1,-Pi,Pi,n),x=-2*Pi..2*Pi); plots:-animate(g,[n],n=1..40);
Here is the animation from the Maple notebook:
Another version
restart; f:=x->piecewise(-Pi<x and x<Pi/2,-1, Pi/2<x and x<1,0,1); assume(n,integer); nmaFourier2:=proc(f,freq,from_,to_,maxN::integer) local n::integer,denomC,denomS,a,b; denomC:=( to_ - from_ ) / 2; denomS:=( to_ - from_ ) / 2; a:=proc(n) int(f(x)*cos(n*freq*x),x=from_..to_) /denomC; end proc; b:=proc(n) int(f(x)*sin(n*freq*x),x=from_..to_) / denomS; end proc; 1/2*a(0) + sum( a(n) * cos(n*freq*x) ,n=1..maxN) + sum( b(n) * sin(n*freq*x) ,n=1..maxN) end proc; plots[setoptions](title=` `, axesfont=[SYMBOL,8] ,font=[COURIER,1], xtickmarks=[seq(evalf(k*Pi/2)=sprintf("%a %s", k/2 ,"pi" ),k= -3..3)], ytickmarks=[-1.0="-1",-0.5="",0.0="0",0.5="",1.0="1"]); B:=array(1..3,1..3); k:=0; for i from 1 to 3 do for j from 1 to 3 do k:=k+1; B[i,j]:=plot({f(x),nmaFourier2(f,1,-Pi,Pi,k)},x=-Pi..Pi,size=[200,100]); end do; end do; plots:-display( B);
restart; v:=1; B:=Matrix(3,3); for i from 1 to 3 do for j from 1 to 3 do v:=v+1; B[i,j]:= plot(x^v,x=-2..2,thickness=3,size=[200,100] ); end do; end do; plots:-display(B);
From book Maple animation by John Putz
plot( sin(x), x=0..2*Pi, xtickmarks=evalf([Pi/2="p/2", Pi="p", 3*Pi/2="3p/2", 2*Pi="2p"]), ytickmarks=[-1,1], axesfont=[SYMBOL,16], labels=["",""] );
From Preben Alsholm
res:=FunctionAdvisor(sin): res2:=op(2,eval(res)): map(print,res2);
or answer by Thomas Richard
> FunctionAdvisor( display, sin );
Use convert(expr,parfrac) or convert(f,fullparfrac)
n := 7; f:=sum('a[k]*b[k]','k'=1..n);
\[ a_{{1}}b_{{1}}+a_{{2}}b_{{2}}+a_{{3}}b_{{3}}+a_{{4}}b_{{4}}+a_{{5}}b_{ {5}}+a_{{6}}b_{{6}}+a_{{7}}b_{{7}} \]
from Serge from the net:
restart; with(geom3d): plane(OYZ,x=0,[x,y,z]): plane(OXZ,y=0,[x,y,z]): plane(OXY,z=0,[x,y,z]): c:=1/2:r:=1/4: L:=combinat[permute]([-c$3,c$3],3): S:=seq(sphere(s||i,[point(A||i,op(op(i,L))),r]),i=1..8): draw([OYZ,OXZ,OXY,S]);
Use evalb(). For example evalb(I*sinh(x)=sin(I*x)); gives true
The above does not always work. Only sure way is to do this
> m1 := exp(I*n*x); m2 := (cos(n*x)+I*sin(n*x)); simplify(m1-m2); simplify(m1-convert(m2,exp));
Function by Robert Israel from the net:
restart; thefacts:= [seq(i!,i=2..20)]: getfacts:= proc(x::{algebraic,series}) local i; if type(x, {`+`,`*`,series}) then map(getfacts,x) elif type(x, fraction) then getfacts(numer(x))/getfacts(denom(x)) elif type(x,`^`) then getfacts(op(1,x))^op(2,x) elif type(x,negint) then -getfacts(-x) elif type(x,posint) then for i from 1 to 19 while irem(x, thefacts[i]) = 0 do od: if i = 1 then x elif thefacts[i-1] = x then ``(i)! else ``(i-1)!*getfacts(x/thefacts[i]) fi else x fi end; getfacts(series(sin(x),x));
\[ \text {series} \left ( x-{\frac {{x}^{3}}{ \left ( \left ( 3 \right ) \right ) !}}+{\frac {{x}^{5}}{ \left ( \left ( 5 \right ) \right ) !}}+O \left ( {x}^{7} \right ) ,x,7 \right ) \]
?updates,maple10
Maple 2020.
restart; PDE := diff(u(x,y), y$2 ) + diff(u(x,y), x$2) = 0; BC:= u(x,0)=0, u(x,100)=100, u(0,y)=0, u(10,y)=0; sol:=pdsolve(PDE,[BC] ,numeric); Error, (in pdsolve/numeric) unable to handle elliptic PDEs
Compare to
restart; PDE := diff(u(x,y), y$2 ) + diff(u(x,y), x$2) = 0; BC:= u(x,0)=0, u(x,100)=100, u(0,y)=0, u(10,y)=0; sol:=pdsolve([PDE,BC]);
\[ u \left ( x,y \right ) =\sum _{n=1}^{\infty }-200\,{\frac { \left ( \left ( -1 \right ) ^{n}-1 \right ) {{\rm e}^{10\,\pi \,n}}\sin \left ( 1/10\,n\pi \,x \right ) \left ( {{\rm e}^{1/10\,n\pi \,y}}- {{\rm e}^{-1/10\,n\pi \,y}} \right ) }{\pi \,n \left ( {{\rm e}^{20\,\pi \,n}}-1 \right ) } } \]
Create a new matrix, by appending some rows of one matrix to rows from another matrix:
restart; with(LinearAlgebra): A:=< <1|2|3> , <4|5|6> >;
\[ \left [ \begin {array}{ccc} 1&2&3\\ \noalign {\medskip }4&5&6 \end {array} \right ] \]
B:=< <7|8|10> , <11|12|13> , <14|15|16> >;
\[ \left [ \begin {array}{ccc} 7&8&10\\ \noalign {\medskip }11&12&13 \\ \noalign {\medskip }14&15&16\end {array} \right ] \]
Now append first row of A to last 2 rows of B
C:=< A[1,1..-1] , B[2..-1,1..-1] >;
\[ \left [ \begin {array}{ccc} 1&2&3\\ \noalign {\medskip }11&12&13 \\ \noalign {\medskip }14&15&16\end {array} \right ] \]
# Now append first column of A to first 2 rows of B A[1..-1,1]; B[1..2,1..-1]; C:=< A[1..-1,1] | B[1..2,1..-1] >;
\[ \left [ \begin {array}{cccc} 1&7&8&10\\ \noalign {\medskip }4&11&12&13 \end {array} \right ] \]
#Now remove the middle row of B B; B:=<B[1,1..-1] , B[-1,1..-1] >;
\[ \left [ \begin {array}{ccc} 7&8&10\\ \noalign {\medskip }14&15&16 \end {array} \right ] \]
#now set the diagonal elements of B to be 0 B:=RandomMatrix(3); for i from 1 to 3 do B[i,i]:=0; end do: B;
\[ B:=\left [ \begin {array}{ccc} 0&99&92\\ \noalign {\medskip }8&0&-31 \\ \noalign {\medskip }69&44&0\end {array} \right ] \] \[ \left [ \begin {array}{ccc} 0&99&92\\ \noalign {\medskip }8&0&-31 \\ \noalign {\medskip }69&44&0\end {array} \right ] \] To find inverse.
restart; with(LinearAlgebra): A:=Matrix( [ [2,0],[4,2] ]); MatrixInverse(A);
\[ \left [ \begin {array}{cc} 1/2&0\\ \noalign {\medskip }-1&1/2 \end {array} \right ] \]
To check that for any matrix A, then A*transpose(A) is always a matrix which is
symmetrical
A:=RandomMatrix(2,3); A.Transpose(A);
\[ A:=\left [ \begin {array}{ccc} 99&44&-31\\ \noalign {\medskip }29&92&67 \end {array} \right ] \]
\[ \left [ \begin {array}{ccc} 99&44&-31\\ \noalign {\medskip }29&92&67 \end {array} \right ] \]
how to create a random lower triangular matrix?
restart; with(LinearAlgebra); A:=RandomMatrix(4,4,outputoptions=[shape=triangular[lower]]);
\[ \left [ \begin {array}{cccc} 67&0&0&0\\ \noalign {\medskip }-31&92&0&0 \\ \noalign {\medskip }44&29&99&0\\ \noalign {\medskip }69&8&27&-4 \end {array} \right ] \]
restart; with(LinearAlgebra); A:=RandomMatrix(5); LinearAlgebra:-Map[(i,j)->evalb(i=j)](x->1,A);
\[ A:= \left [ \begin {array}{ccccc} 1&-98&-76&-4&29\\ \noalign {\medskip }-38& 1&-72&27&44\\ \noalign {\medskip }-18&57&1&8&92\\ \noalign {\medskip }87& 27&-32&1&-31\\ \noalign {\medskip }33&-93&-74&99&1\end {array} \right ] \] \[ \left [ \begin {array}{ccccc} 1&-98&-76&-4&29\\ \noalign {\medskip }-38& 1&-72&27&44\\ \noalign {\medskip }-18&57&1&8&92\\ \noalign {\medskip }87& 27&-32&1&-31\\ \noalign {\medskip }33&-93&-74&99&1\end {array} \right ] \]
eq:=3*x^3+2*x^2+x+5=0; s:=[evalf(solve(eq,x))]; mul(s[i],i=1..nops(s));
Gives
restart; eq:=3*x+4*y+2*z=10; plot3d(solve(eq,z),x=-5..5,y=-5..5,axes=normal);
One can also use impliticplot3d
restart; with(plots): implicitplot3d(3*x+4*y+2*z=10, x=-5..5,y=-5..5, z=-20..20,axes=normal);
From http://www.mapleprimes.com/questions/40470-Trigonometric-Function-To-Sinc-Function
Maple doesn’t have a sinc function. If you mean the function sinc(x) = sin(x)/x, you could say something like
> eval(expr, {sin = (x -> x*sinc(x)), cos = (x -> (x+Pi/2)*sinc(x+Pi/2)), tan = (x -> x*sinc(x)/(x+Pi/2)/sinc(x+Pi/2))});
restart; with(LinearAlgebra): A:=Matrix([[1,0,1,0,1],[0,1,0,1,0]]); NullSpace(A); ColumnSpace(A);
Go to tools->optiopn, and Display, and select Maple notation for input display.
solve(x^2-sin(x),x); RootOf(-sin(_Z)+_Z^2) allvalues(%); RootOf(-sin(_Z)+_Z^2, 0.), RootOf(-sin(_Z)+_Z^2, .8767262154) evalf(%); 0., .8767262154
Use Map with filter
A:=< 1,2,3;4,5,6;7,8,9>; LinearAlgebra:-Map[(i,j)->evalb(i=j)](x->x+1,A);
Go to http://www.maplesoft.com/support/help/search.aspx
and type say updates,Maple17,DE in the small box there.
From http://www.mapleprimes.com/questions/201092-How-To-Insert-New-Paragraph-On-Its-Own by Carl Love:
I use these special keystrokes constantly in my Maple worksheet typing: Ctrl-J: Insert execution group below cursor. Ctrl-K: Insert execution group above cursor. Ctrl-T: Switch from executable code mode to text mode (for entering extended formatted comments). Ctrl-M: Switch from text mode to executable code mode. Shift-Enter (or Shift-Return): Begin a new line in the same execution group. Func-3: Split execution group into two (at cursor). Func-4: Join cursor execution group with execution group below.
Use the read command, as in read "mycode.mpl" where mycode.mpl is plain text file that
contains maple code
New packages are module, which allows using packageName:-function() since it is easier.
Old packages use tables which needs packageName[function]() which is not
common.
To find if package is based on module or not, use the command
type(combstruct,'`module`');
This will return true or false. To know if name is package use the command
type(combstruct,'package');
file_name :=StringTools:-SubstituteAll(file_name,":","-");
restart; c:= i->([i/(1+i),0],1/(1+i)): d:= i->([1,1/i],1/i): geometry:-circle(c1,[geometry:-point(o,2/3,0),1/3],[x,y]): geometry:-circle(c2,[geometry:-point(o,1,1),1],[x,y]): geometry:-intersection(o,c1,c2,[u,v]): plots:-display(plottools:-circle(c(2)),plottools:-circle(d(1)),geometry:-draw(o));
To know more about the intersection, use this:
geometry:-detail(o);
Use symbolic option
restart; simplify(ln(3^x/2^y) =ln(n),symbolic);
How to convert \[ \frac {3+2\sinh (x)^2}{\sinh (x)^2\tanh (x)} \] to \[ 3 \coth ^3(x)-\coth (x) \]
restart; e := (3+2*sinh(x)^2)/(sinh(x)^2*tanh(x)); expand(student[changevar](sinh(x)^2=tanh(x)^2/(1-tanh(x)^2),e));
restart; try fd :=-1; fd := fopen("C:\\output3.txt",APPEND,TEXT); catch: print(`Unable to open file, error is`); print(StringTools:-FormatMessage(lastexception[2])); end try: if not(evalb(fd=-1)) then #file open ok str:="hello world"; try fprintf(fd,"%s\n",str); catch: print(`failed to append to file, error is`); print(StringTools:-FormatMessage(lastexception[2])); finally: close(fd); end try; fi:
To find in which library a command is do
with(LibraryTools); FindLibrary('int',all); #find which library command int is in "C:\Program Files\Maple 18\lib\update.mla", "C:\Program Files\Maple 18\lib\DEsAndMathematicalFunctions18.mla", "C:\Program Files\Maple 18\lib\maple.mla"
To get content of library do
restart; with(LibraryTools): LibLocation:=cat(kernelopts(mapledir),"/lib/maple.mla"); c:=ShowContents(LibLocation);
Then can use this to print the name of each symbol/command, and then use whattype command to find its type
seq(c[i,1],i=1..20);
To get list of Maple kernel builtin commands and symbols, use this. Written by Acer from Maple prime site:
restart: interface(warnlevel=0): started := false: T := 'T': for i from 1 to 1000 do f := eval(parse(cat("proc() option builtin=",i,"; end proc"))); p := (s->StringTools:-Take(s,StringTools:-Search(";",s)-1))(convert(eval(f),string)[26..]); if not type(parse(p),posint) then T[i] := p; started := true; else if started then i:=1000; next; end if; end if; end do: i; [ entries(T,nolist) ]; nops(%);
The above gives on Maple 18.02 the following
["crinterp", "equation", "`{}`", "even", "debugopts", "embedded_imaginary", "define_external", "embedded_real", "coeff", "cx_zero", "coeffs", "embedded_axis", "conjugate", "constant", "convert", "cx_infinity", "dlclose", "identical", "divide", "hfloat", "`done`", "function", "`$`", "fraction", "denom", "float", "degree", "finite", "disassemble", "extended_rational", "diff", "extended_numeric", "frem", "`union`", "frontend", "upperbound", "exports", "writeto", "factorial", "`xor`", "evalgf1", "type", "expand", "typematch", "entries", "unames", "evalb", "unbind", "`evalf/hypergeom/kernel`", "atomic", "hfarray", "anything", "hastype", "complex", "has", "boolean", "goto", "`:-`", "gmp_isprime", "`!`", "genpoly", "anyfunc", "gc", "algebraic", "SFloatMantissa", "ssystem", "Scale10", "`stop`", "Scale2", "sort", "SearchText", "`[]`", "`~`", "`subset`", "~Array", "subsindets", "~Matrix", "streamcall", "~Vector", "subs", "Unordered", "table", "ToInert", "system", "_hackwareToPointer", "substring", "UpdateSource", "subsop", "_maplet", "trunc", "_jvm", "`kernel/transpose`", "_treeMatch", "tcoeff", "_savelib", "taylor", "abs", "rtable_num_dims", "addressof", "rtable_num_elems", "_unify", "rtable_options", "_xml", "rtable_redim", "`and`", "rtable_scale", "andmap", "rtable_scanblock", "alias", "rtable_size", "anames", "rtable_sort_indices", "assign", "savelib", "assemble", "rtable_zip", "array", "select", "appendto", "searchtext", "cat", "series", "callback", "selectremove", "bind", "sign", "attributes", "setattribute", "ormap", "ArrayOptions", "order", "Array", "parse", "`**`", "overload", "`*`", "`::`", "numer", "CopySign", "numelems", "`^`", "`or`", "`||`", "op", "nops", "seq", "normal", "time", "`not`", "piecewise", "numboccur", "`?[]`", "userinfo", "modp2", "inner", "mods", "timelimit", "mvMultiply", "traperror", "negate", "rtable_normalize_index", "call_external", "rtable_is_zero", "assigned", "rtable_indfns", "evalf", "rtable_histogram", "eval", "evaln", "rtable_eval", "truefalse", "evalhf", "rtable_convolution", "tabular", "mul", "rtableInfo", "zppoly", "`if`", "rtable", "uneval", "remove", "sfloat", "rhs", "specfunc", "readlib", "string", "reduce_opr", "symbol", "ASSERT", "`?()`", "realcons", "TRACE", "`quit`", "relation", "_local", "pointto", "sequential", "add", "print", "set", "SFloatExponent", "iolib", "radical", "SDMPolynom", "`int/series`", "protected", "Record", "irem", "procedure", "Re", "iquo", "poszero", "isqrt", "real_infinity", "RETURN", "is_gmp", "ratpoly", "`+`", "lcoeff", "rational", "OrderedNE", "kernelopts", "range", "Object", "NumericEventHandler", "icontent", "numeric", "NumericStatus", "igcd", "odd", "NumericClass", "ilog10", "nonpositive", "NumericEvent", "ilog2", "nonreal", "`implies`", "posint", "NameSpace", "indets", "positive", "NextAfter", "indices", "polynom", "MPFloat", "`intersect`", "pos_infinity", "MorrBrilCull", "`<`", "member", "neg_infinity", "Im", "maxnorm", "name", "`<>`", "max", "negint", "`<=`", "map2", "negative", "modp1", "nonnegative", "FromInert", "modp", "negzero", "EqualStructure", "`minus`", "nonposint", "`>=`", "min", "nonnegint", "`>`", "DefaultUnderflow", "lexorder", "imaginary", "`=`", "lhs", "indexable", "ERROR", "ldegree", "indexed", "EqualEntries", "length", "integer", "macro", "list", "DEBUG", "map", "literal", "`..`", "lowerbound", "`module`", "Default0", "lprint", "moduledefinition", "DefaultOverflow"] 296
This one has one solution
eq:=diff(u(z),z$2)+(k-1)*diff(u(z),z)/z+lambda*exp(u(z))=0; sol:=dsolve({subs({k=1,lambda=2},eq),u(0)=1,u(1)=0},numeric,u(z), method=bvp[midrich],'abserr'=0.001); plots[odeplot](sol);
This solved coupled ODE’s, so there are 2 solutions. Say \(x_1(t)\) and \(x_2(r)\), It is a little tricky to plot all solutions generated, but here is an example
restart; R := 0.4; px := 32000; Mm := 0.1; Ds := 9; DO2 := 7.2; YXS := 0.3; KS := 10; Sp := 30; Cb := 8; KO2 := 0.2; R0 := 0.000001; YXO := 0.42857; Vs := px*1/YXS*(Mm*x2(r))/(KS + x2(r))*x1(r)/(KO2 + x1(r)); Vo := px*1/YXO*(Mm*x2(r))/(KS + x2(r))*x1(r)/(KO2 + x1(r)); eqs := diff(x1(r),r$2) + 2/r*diff(x1(r),r)= Vo/DO2, diff(x2(r),r$2) + 2/r* diff(x2(r),r)= Vs/Ds; ic:=D(x1)(R0)=0,x1(R) = Cb,D(x2)(R0)= 0, x2(R) = Sp; sol:=dsolve({eqs,ic},numeric,{x1(r),x2(r)},'abserr'=.52,'maxmesh'=1000,output=listprocedure);
And now to plot do
x1Sol:=rhs(sol[2]); plot(x1Sol(r),r=0..0.4); x2Sol:=rhs(sol[4]); plot(x2Sol(r),r=0..0.4);
This below by Axel Vogt posted on sci.math.symbolic which does a nice job of formatting
output to specific width.
split_for_print:=proc(expr, len) # expr = some Maple expression # len = length to split with line breaks local L,s,tmp,j; s:=convert(expr, string); L:=[StringTools:-LengthSplit(s, len)]; for j from 1 to nops(L) do # if j = nops(L) then printf("%s ;", L[-1]) if j = nops(L) then printf("%s", L[-1]) else printf("%s\\\n", L[j]); end if; end do: end proc; evalf[100](Pi); split_for_print(%, 40); 3.14159265358979323846264338327950288419\ 7169399375105820974944592307816406286208\ 998628034825342117068
for VIM
in vim, type set syntax=maple after putting the file maple.vim in ~/.vim/syntax/maple.vim.
I found maple.vim in above link.
For Maple IDE
use packages(); to find what packages loaded. use unwith to remove package
packages(); [] with(DynamicSystems): packages(); [DynamicSystems] unwith(DynamicSystems); packages(); []
restart in separate execution group
with inside proc(). Use uses instead.
assume( A::AndProp(NonZero,constant) );
Now can use is(A,constant);
check for ‘=‘ as follows
eq:= x=1; whattype(eq); # `=` if whattype(eq) = `=` then print("yes"); else print("no"); fi; "yes"
check for ‘set‘ as follows
eq:= {diff(y(x),x)=1,x(0)=1}; if whattype(eq) = `set` then print("yes"); else print("no"); fi; "yes"
I could only find a way to export to eps
plotsetup(default): plotsetup(postscript, plotoutput=`t.eps`, plotoptions=`color,portrait,height=300`); plot(sin(x),x=-Pi..Pi,'gridlines'); plotsetup(default):
Make sure not to put : at the end of the plot command! else it will not be exported. It has
to end with ;
This will same it to t.eps in the currentdir() location. Then used ps2pdf t.eps t.pdf
to convert it to PDF. Or just ps2pdf t.eps it will automatically create t.pdf
Or ps2pdf -dCompatibilityLevel=1.4 t.eps but may it is best to do
ps2pdf -dCompatibilityLevel=1.4 -dEmbedAllFonts=true t.eps
Also try adding
-dPDFSETTINGS=/printer
to the above. This tells it to optimize it for printing.
Another example of a direction field for an ODE
plotsetup(postscript, plotoutput=`t0.eps`, plotoptions=`color,portrait, height=300` ); ode:= diff(y(x),x) = 3*x^2 - 1; DEtools:-DEplot( ode, y(x), x=-2..2, [y(0) = 0], y=-2..2, linecolour=red, color = blue, stepsize=.05,arrows=MEDIUM ); plotsetup(default);
To find roots of \( (3+4 i)^{1/3}\), do
fsolve(z^3=(3+4*I),z); #gives -1.26495290635775+1.15061369838445*I, -.363984239564424-1.67078820068900*I, 1.62893714592218+.520174502304545*I
A:= Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 2}); f:=x->`if`(x<>0,x*LinearAlgebra:-IdentityMatrix(2),0*Matrix(2)); B:=map(f,A);
Which gives
\[ \left [ \begin {array}{cc} \left [ \begin {array}{cc} 0&0 \\ \noalign {\medskip }0&0\end {array} \right ] & \left [ \begin {array} {cc} 0&0\\ \noalign {\medskip }0&0\end {array} \right ] \\ \noalign {\medskip } \left [ \begin {array}{cc} 0&0 \\ \noalign {\medskip }0&0\end {array} \right ] & \left [ \begin {array} {cc} 2&0\\ \noalign {\medskip }0&2\end {array} \right ] \end {array} \right ] \]
now
r:=Matrix(convert(B,listlist))
Gives
\[ \left [ \begin {array}{cccc} 0&0&0&0\\ \noalign {\medskip }0&0&0&0 \\ \noalign {\medskip }0&0&2&0\\ \noalign {\medskip }0&0&0&2\end {array} \right ] \]
Maple has a simple but easy to use pattern matching, which works well. Here are some example. For each case, will show what pattern to detect and how to do it. I am still not very good at pattern matching in Maple and will need to make improvement in this with time.
Detect \(\sqrt {x y}\) in expression.
restart; expr:= sin(x)*sqrt(x*y); if patmatch(expr,a::anything*(b::anything*x*y)^(c::anything),'la') then assign(la); if c =1/2 or c=-1/2 then print("found sqrt(x*y)"); else print("did not find sqrt(x*y)"); fi; fi;
But if the expression was \(\sin (x)\sqrt {x y}+3\) then the above would fail, because there are a term after \(\sqrt {x y}\), so the pattern has to change to
restart; expr:= sin(x)*sqrt(x*y)+3; if patmatch(expr,a::anything*(b::anything*x*y)^(c::anything)+d::anything,'la') then assign(la); if c =1/2 or c=-1/2 then print("found sqrt(x*y)"); else print("did not find sqrt(x*y)"); fi; fi;
There was a case where I wanted to detect form \(f(x) g(\frac {y}{x})\), i.e. \(f(x)\) can be any expression which is function of \(x\) only (it can be constant also) multiplied by a function whose argument must be \(\frac {y}{x}\) or a constant multiplied by \(\frac {y}{x}\).
This means something like \(x g(\frac {y}{x})\) or \(x^2 e^{3 \frac {y}{x}}\) or \(f(x) \sin {\frac {y}{x}}\) or \(\cos {\frac {y}{x}}\) where in this last case \(f(x)=1\) which is allowed.
TO FINISH.
use trigsubs, very useful command. For example
trigsubs(cos(theta)^3)
Gives
\[ [1/2\,\cos \left ( \theta \right ) +1/2\,\cos \left ( 2\,\theta \right ) \cos \left ( \theta \right ) ,1/4\,\cos \left ( 3\,\theta \right ) +3/4\,\cos \left ( \theta \right ) ] \]
Given \(f(x,y,z)=x^2 z+y^3 z^2-xyz\) we want to find its directional derivative along the vector \(n\).
One way
n:=<-1,0,3>; g:=VectorCalculus[Gradient](x^2*z+y^3*z^2-x*y*z, [x,y,z]); Student[VectorCalculus][DotProduct](g,n/LinearAlgebra[Norm](n,2))
Gives
\[ -{\frac { \left ( 2\,xz-yz \right ) \sqrt {10}}{10}}+{\frac { \left ( 6\,{y}^{3}z+3\,{x}^{2}-3\,xy \right ) \sqrt {10}}{10}} \]
Another is
Student[MultivariateCalculus][DirectionalDerivative](x^2*z+y^3*z^2-x*y*z, [x,y,z], [-1,0,3]);
Gives the same result.
For simple variable, use assigned
restart; x:=10: assigned(x) true assigned(y) false
For a field in table do
restart; A:=table(["x"=10,"y"=20]): assigned(A["x"]) true assigned(A["z"]) false
For field in Record, I do not know how yet, other than using try catch, as assigned does not seem to work for Record fields.
restart; A:=Record('x'=10,'y'=20); try assigned(A:-x) catch: print("no such field in record") end try; true try assigned(A:-z) catch: print("no such field in record") end try; "no such field in record"
given \[ {\mathrm e}^{\frac {2 \ln \left (\sqrt {p^{2}+1}+p \right )+2 \ln \left (a \right )+\ln \left (p^{2}+1\right ) a}{2 a}}+{\mathrm e}^{3 x} \]
simplify(expr) does not work. So tried subsindets
restart; expr := exp((2*ln(sqrt(p^2 + 1) + p) + 2*ln(a) + ln(p^2 + 1)*a)/(2*a))+ exp(3*x); subsindets(expr,'specfunc( anything, exp )',f->(`if`(has(op(1,f),'ln'),expand(f),f)))
\[ \left (\sqrt {p^{2}+1}+p \right )^{\frac {1}{a}} a^{\frac {1}{a}} \sqrt {p^{2}+1}+{\mathrm e}^{3 x} \]
It is possible to also try simplify(expr,exp) in some cases, but for the above example, this
did not work, i.e. it did not simplify it.
Update december 2023. Trying Maple 2023.2.1, it simplifies the above using
simplify(expr,exp)
restart; expr := exp((2*ln(sqrt(p^2 + 1) + p) + 2*ln(a) + ln(p^2 + 1)*a)/(2*a))+ exp(3*x); simplify(expr,exp)
\[ sqrt{p^{2}+1}\, \left (\sqrt {p^{2}+1}+p \right )^{\frac {1}{a}} a^{\frac {1}{a}}+{\mathrm e}^{3 x} \]
And
restart; expr:=exp(ln(x)+ln(y)); simplify(expr)
\[ x y \]
Given \[ \left [\begin {array}{cccc}1 & -1 & 0 & 2 \\1 & 2 & 2 & -2 \\0 & 2 & 3 & -1 \end {array}\right ] \]
Find its Null, Row and Column space basis vectors.
restart; A:=Matrix([[1,-1,0,2],[1,2,2,-2],[0,2,3,-1]]); LinearAlgebra:-NullSpace(A)
\[ \left \{ \left [\begin {array}{c}0 \\2 \\-1 \\1 \end {array}\right ] \right \} \]
restart; A:=Matrix([[1,-1,0,2],[1,2,2,-2],[0,2,3,-1]]); LinearAlgebra:-RowSpace(A)
\[ \left [\left [\begin {array}{cccc}1 & 0 & 0 & 0 \end {array}\right ], \left [\begin {array}{cccc}0 & 1 & 0 & -2 \end {array}\right ], \left [\begin {array}{cccc}0 & 0 & 1 & 1 \end {array}\right ]\right ] \]
restart; A:=Matrix([[1,-1,0,2],[1,2,2,-2],[0,2,3,-1]]); LinearAlgebra:-ColumnSpace(A)
\[ \left [\left [\begin {array}{c}1 \\0 \\0 \end {array}\right ], \left [\begin {array}{c}0 \\1 \\0 \end {array}\right ], \left [\begin {array}{c}0 \\0 \\1 \end {array}\right ]\right ] \]
Given \[ \left [\begin {array}{cccc}1 & -4 & -3 & -7 \\2 & -1 & 1 & 7 \\1 & 2 & 3 & 11 \end {array}\right ] \]
Find the new form after Gaussian elimination
restart; A:=Matrix([[1,-4,-3,-7],[2,-1,1,7],[1,2,3,11]]); LinearAlgebra:-GaussianElimination(A);
\[ \left [\begin {array}{cccc}1 & -4 & -3 & -7 \\2 & -1 & 1 & 7 \\1 & 2 & 3 & 11 \end {array}\right ] \]
Given matrix \[ \left [\begin {array}{ccc}5 & 2 & 18 \\0 & 1 & 4 \\4 & 1 & 12 \end {array}\right ] \]
Find its Reduced Echelon form.
restart; A:=Matrix([[5,2,18],[0,1,4],[4,1,12]]); Student:-LinearAlgebra:-ReducedRowEchelonForm(A)
\[ \left [\begin {array}{ccc}1 & 0 & 2 \\0 & 1 & 4 \\0 & 0 & 0 \end {array}\right ] \]
Another option is
restart; A:=Matrix([[5,2,18],[0,1,4],[4,1,12]]); MTM:-rref(A)
\[ \left [\begin {array}{ccc}1 & 0 & 2 \\0 & 1 & 4 \\0 & 0 & 0 \end {array}\right ] \]
Given matrix \[ \left [\begin {array}{cc}1 & 1 \\2 & 3 \\4 & 5 \end {array}\right ] \] How to add row \[ [a, b] \] to end of the matrix?
restart; A:=Matrix([[1,1],[2,3],[4,5]]); the_row:=convert([a,b],Vector['row']); ArrayTools:-Concatenate(1,A,the_row);
\[ \left [\begin {array}{cc}1 & 1 \\2 & 3 \\4 & 5 \\a & b \end {array}\right ] \]
Use LinearAlgebra:-Adjoint and then transpose the result. Since the Adjoint is the
transpose of the cofactor.
Given
\[ \left [\begin {array}{ccc} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 10 \end {array}\right ] \]
then
restart; A:=Matrix([[1,2,3],[4,5,6],[7,8,10]]); LinearAlgebra:-Transpose(LinearAlgebra:-Adjoint(A))
\[ \left [\begin {array}{ccc} 2 & 2 & -3 \\ 4 & -11 & 6 \\ -3 & 6 & -3 \end {array}\right ] \]
When finding eigenvectors of matrix, using LinearAlgebra, the vectors are not normalized.
How to normalized them so the length is one?
One way is
restart; LA:=LinearAlgebra; Sx:=Matrix([[0,1,0],[1,0,1],[0,1,0]]); #this finds eigenvectors in v lam,v:=LA:-Eigenvectors(Sx); #this normalize it B:=map(n -> v[.., n]/norm(v[.., n], 2), [$1..LA:-RowDimension(v)]): B:=`<|>`(op(B)); #this converts the list back to matrix.
\[ v=\left [\begin {array}{ccc} -1 & 1 & 1 \\ 0 & \sqrt {2} & -\sqrt {2} \\ 1 & 1 & 1 \end {array}\right ] \]
\[ B=\left [\begin {array}{ccc} -\frac {\sqrt {2}}{2} & \frac {1}{2} & \frac {1}{2} \\ 0 & \frac {\sqrt {2}}{2} & -\frac {\sqrt {2}}{2} \\ \frac {\sqrt {2}}{2} & \frac {1}{2} & \frac {1}{2} \end {array}\right ] \]
expr:=`ℏ`*x
gives
\[ \hslash x \]
Notice, the ; is needed. This `&hbar`*x will not work. It must be `ℏ`*x
First example
restart; VectorCalculus:-SetCoordinates( 'cartesian'[x,y,z] ); F:=VectorCalculus:-VectorField(<y,-x,0>);
\begin {align*} F &= y \mathbf {\bar {e}_{x}}-x \mathbf {\bar {e}_{y}} \end {align*}
And now
VectorCalculus:-Curl(F);
\[ -2 \mathbf {\bar {e}_{z}} \]
Second example
restart; VectorCalculus:-SetCoordinates( 'cartesian'[x,y,z] ); F:=VectorCalculus:-VectorField(<y*z^2,x*z^2+2,2*x*y*z-1>);
\begin {align*} F &= y \,z^{2} \mathbf {\bar {e}_{x}}+(x \,z^{2}+2) \mathbf {\bar {e}_{y}}+(2 x y z -1) \mathbf {\bar {e}_{z}} \end {align*}
And now
VectorCalculus:-Curl(F);
\[ 0 \] Since Curl is zero, field is conservative.
Third example, in cylinderical coodinates
restart; VectorCalculus:-SetCoordinates( 'cylindrical'[rho,phi,z] ); F:=VectorCalculus:-VectorField(<0,-rho,2>);
\begin {align*} F &= -\rho \mathbf {\bar {e}_{\phi }}+2 \mathbf {\bar {e}_{z}} \end {align*}
And now
VectorCalculus:-Curl(F);
\[ 2 \mathbf {\bar {e}_{z}} \]
Use Student:-LinearAlgebra:-GaussJordanEliminationTutor( A, output=steps )
Where \(A\) is your augmented matrix.
Do not use the Maple command LinearAlgebra:-ColumnSpace for this. it gives the
columns in the RREF. The correct way is to obtain the corresponding columns of
the pivot columns in the original matrix \(A\). Hence use the command Basis like
this
A:=Matrix([[1,0,0],[1,1,1]]); LinearAlgebra:-Basis([seq(A[..,i],i=1..LinearAlgebra:-ColumnDimension(A) )]);
Which gives
\[ \left [\left [\begin {array}{c} 1 \\ 1 \end {array}\right ], \left [\begin {array}{c} 0 \\ 1 \end {array}\right ]\right ] \]
If you use ColumnSpace command you’ll get this
A:=Matrix([[1,0,0],[1,1,1]]); LinearAlgebra:-ColumnSpace(A);
\[ \left [\left [\begin {array}{c} 1 \\ 0 \end {array}\right ], \left [\begin {array}{c} 0 \\ 1 \end {array}\right ]\right ] \]
These are different. Basis is the correct command to use, which matches the standard
definition in textbooks.
For integration do
Student:-Calculus1:-ShowSolution(Int(x*sin(x),x));
The steps are displayed. This does not work all the time. For example
integrand:=x*y(x)*diff(y(x),x$2)+x*(diff(y(x),x))^2-y(x)*diff(y(x),x); Student:-Calculus1:-ShowSolution(Int(integrand,x));
gives
Error, (in Student:-Calculus1:-ShowSolution) unable to determine which calculus operation is being applied in this problem; you can provide this information as the 2nd argument on your call to Rule or Hint
For differential equations, support is limited but these are the steps
restart; ode:=diff(y(x),x)=sin(x); Student:-ODEs:-ODESteps(ode)
Prints the steps. If IC is there, then
restart; ode:=diff(y(x),x)=sin(x); ic:=y(0)=1; Student:-ODEs:-ODESteps([ode,ic])
Use FileTools:-ListDirectory
dir_name:="C:/tmp"; currentdir(dir_name); #cd to directory files_to_process := FileTools:-ListDirectory(dir_name,'all','returnonly'="*.tex"): numelems(files_to_process) 100
In the above, files_to_process is a list of the files in the current folder with extension
.tex
There was a case when I needed to delete lines from text file that contains a say "foo" as an example.
This is what I did. use readline to read the lines, check, and if the line contains "foo" skip, else write the line to a temporary file. At the line, use Rename to rename the temporary file to the file being read.
dir_name:="C:/tmp"; currentdir(dir_name); tmp_file_name := "TMP.txt"; source_file_name := "source.txt"; file_id := fopen(tmp_file_name,WRITE): line := readline(source_file_name): while line<>0 do if not StringTools:-Has(line,"foo") then fprintf(file_id,"%s\n",line); fi; line := readline(source_file_name): od: fclose(file_id); FileTools:-Rename(tmp_file_name,source_file_name,force=true);
Given \(9 x^{5}+4 x^{4}+3 x^{3}+x^{2}+x +1\) how to truncate it, so that all terms of \(x^3\) and higher are removed?
This can be done as follows
restart; p:=1+x+x^2+3*x^3+4*x^4+9*x^5; simplify(p,{x^3=0})
\[ x^{2}+x +1 \]
Sometimes it is useful to make a small local piece of code inside a proc, with its
own local variables that do not interfer with the variables of the proc. In Ada,
this is done using declare clause for example. In Maple on can do the same as
follows
restart; foo:=proc() local n; n:=10; proc() local n; n:=99; print("inside inner proc, n=",n); end proc(); print("n=",n); end proc; foo();
Which prints
"inside inner proc, n=", 99 "n=", 10
Notice the end of the inner anonymous proc above. It has end proc(); and not end proc;
as normal proc. This defines the inner proc and calls it at same time. All the local variables
inside the anonymous proc only exist inside that proc and go away after the call, and they
do not interfer with the outer proc variables. This way we can declare temporary variables
and use them where they are needed.
There was a case where I was making lots of calls from many places to pne specific proc inside a module. I did not want to keep using the long name each time.
the command alias did not work. After some trial and error, found that using use works.
Here is the solution. First this is the original layout
restart; A:=module() export B:=module() export foo:=proc() print("in A:-B:-foo()"); end proc; end module; export C:=module() export boo:=proc() print("in A:-C:-boo()"); A:-B:-foo(); end proc; end module; end module;
In the above, the goal is replace A:-B:-foo(); with just foo() and have it bind to
A:-B:-foo(); automatically.
This is done by modifying the above to
restart; A:=module() export B:=module() export foo:=proc() print("in A:-B:-foo()"); end proc; end module; use foo=A:-B:-foo in #add this line here export C:=module() export boo:=proc() print("in A:-C:-boo()"); foo(); #now can just use the short name end proc; end module; end use; #add this line here. end module;
Wrapping the whole module where the short name is used worked.
Any module that needs to use the short name, can do the same. This solved the problem.
I had case where there is list of Objects, and wanted to removed duplicate entries in the list based on if some field is the same among the objects.
This can be done using the command ListTools:-MakeUnique and using a proc which
checks for the condition. In this example, we want to remove objects where the field age in
each object is the same.
restart; #create data type module person_type() option object; export name::string:="me"; export age:=25; end module: #make few objects o1:=Object(person_type); o2:=Object(person_type); o3:=Object(person_type); o3:-age:=46; o4:=Object(person_type); #make list of them list_of_people:=[o1,o2,o3,o4]; nops(list_of_people); #this will print 4 #now delete object if age is same list_of_people:=ListTools:-MakeUnique( list_of_people, 1, proc(a,b) evalb(a:-age=b:-age); end proc ); nops(list_of_people); #this will print 2
Converting a list of Vectors to set will not remove duplicates, as each Vector occupies
different memory address, even if the structure is the same. To remove duplicate vector, use
ListTools:-MakeUnique as follows
restart; my_list:=[Vector([1,0]),Vector([1,0]),Vector([2,0])]; convert(my_list,set); #this will still show the 3 vectors. ListTools:-MakeUnique(my_list,1,proc(a,b) LinearAlgebra:-Equal(a,b) end proc) #now only 2 vectors will remain. Duplicate one was removed
Gives a rational function in \(x\), such as \[ \frac {1}{10 \left (x -4\right ) \left (x -5\right )^{3}} \] How to find all its poles which are \(x=4\) and \(x=5\) and the order of each pole which will be \(1\) and \(3\) in this example?
Using sqrfree as follows
restart; get_poles_and_order:=proc(r_in,x::symbol)::list; local r:=r_in,N::posint; local the_poles::list; local item; r:=normal(r); if not type(r,'ratpoly'(anything,x)) then error("Not be a polynomial or a rational function in ",x) fi; the_poles := sqrfree(denom(r),x); the_poles := the_poles[2,..]; #we do not need the overall factor for N,item in the_poles do the_poles[N]:=[solve(item[1]=0,x),item[2]]; od; return the_poles; end proc:
The above proc get_poles_and_order returns back a list of lists. Each sublist has its first
entry the pole and the second entry the order.
Here are some examples
r:=1/(10*(x-4)*(x-5)^3); get_poles_and_order(r,x) #[[4, 1], [5, 3]]
The above says there is a pole at \(x=4\) of order 1 and pole at \(x=5\) of order 3.
Doing series in Maple with specific order value, the number of terms generated ofcourse
depends on the function. I had need to have the series generated always with same number
of terms. I could not find an option in Maple to do that. This function does this. It keeps
finding the series for the function with increasing order until the number terms that comes
out is what requested. There is an upper limit that can be changed if needed to protect
against pathological cases.
restart; get_series_by_terms:=proc(expr,x::symbol,at::numeric,number_terms_needed::posint) local keep_running::boolean:=true; local current_order::integer:=0; local MAX_ORDER_TO_TRY::posint:=100; #change as needed local result; do current_order := current_order+1; result := convert(series(expr,x=at,current_order),polynom); if nops(result) >= number_terms_needed or current_order>MAX_ORDER_TO_TRY then keep_running:=false; fi; until keep_running=false; return result; end proc:
And now
get_series_by_terms(sin(x),x,0,10)
returns \begin {equation} \begin {aligned} & x -\frac {1}{6} x^{3}+\frac {1}{120} x^{5}-\frac {1}{5040} x^{7}+ \frac {1}{362880} x^{9} -\\ & \frac {1}{39916800} x^{11}+\frac {1}{6227020800} x^{13}-\frac {1}{1307674368000} x^{15}+\\ & \frac {1}{355687428096000} x^{17}-\frac {1}{121645100408832000} x^{19} \end {aligned} \end {equation}
Given a parent module \(A\) and inside it there are two child modules (local modules) with names
say \(B\) and \(C\). To call a proc foo inside \(B\) from another proc inside \(C\), the proc foo has to be
exported. But the module \(B\) does not have to be exported, if we make sure to use B:-foo()
call instead of full name A:-B:-foo() call.
So make sure to use child:-proc() from other sibilings to avoid having to make each child
exported. Making children exported means they can be seen and called directly from outside
the parent which his not what we want.
Here is an example
restart; A:=module() #parent export main:=proc() C:-foo(); end proc; local B:=module() #child export foo:=proc() print("in A:-B:-foo() proc"); end proc; end module; local C:=module() #child export foo:=proc() print("in A:-C:-foo(). About to call A:-B:-foo()"); B:-foo(); #do this and NOT A:-B:-foo() end proc; end module; end module;
and now
A:-main() "in A:-C:-foo(). About to call A:-B:-foo()" "in A:-B:-foo() proc"
If instead we have written A:-B:-foo() in the above call, then Maple will complain with
the error Error, (in foo) module does not export `B`
Maple’s command Value(Time()) returns 13 digits number, which is number of milliseconds
from epoch. I wanted this value to be in seconds, to match the file changed time
from FileTools[Status]("A.txt" ) which uses seconds and not milliseconds. I
could not find an option to tell Date or Time to do this. Here is one way to do
this.
r:=Value(Time()); #r := 1652677498870 length(r); #13 r:=convert(r, base, 10); r:=ListTools:-Reverse(r); r:=r[1..-4]; #remove last 3 digits nops(r); r:=parse(cat(op(r))) #r := 1652677498 length(r); #10
This can be made into a function
get_time_in_seconds:=proc()::integer; local r; r:=Value(Time()); r:=convert(r, base, 10); r:=ListTools:-Reverse(r); r:=r[1..-4]; r:=parse(cat(op(r))); return r; end proc; get_time_in_seconds() #1652679222
I noticed that Maple returns the summation index variable using leading underscore as in _n
or _m which makes the latex looks not as good. Here is an example
restart; dsolve(diff(y(x),x$2)+diff(y(x),x)+y(x)=0,y(x),'formal_series'); y(x) = _C1*Sum((-1/2 - sqrt(3)*I/2)^_n*x^_n/_n!, _n = 0 .. infinity) + _C2*Sum((-1/2 + sqrt(3)*I/2)^_n*x^_n/_n!, _n = 0 .. infinity)
The latex of the above is \[ y \left (x \right ) = \textit {\_C1} \left (\sum _{\textit {\_n} =0}^{\infty }\frac {\left (-\frac {1}{2}-\frac {i \sqrt {3}}{2}\right )^{\textit {\_n}} x^{\textit {\_n}}}{\textit {\_n} !}\right )+\textit {\_C2} \left (\sum _{\textit {\_n} =0}^{\infty }\frac {\left (-\frac {1}{2}+\frac {i \sqrt {3}}{2}\right )^{\textit {\_n}} x^{\textit {\_n}}}{\textit {\_n} !}\right ) \]
Not seeing an option to change _n to n, I wrote the following function which takes in the
solution, use subsindets and remove the leading underscore.
This is the above example showing how to use the function
restart; fix_summation_index:=proc(expr) local fix_it:=proc(the_sum) local the_letter::symbol,the_new_letter::symbol,the_letter_as_string::string; the_letter:= op([2,1],the_sum); the_letter_as_string:=String(the_letter); if the_letter_as_string[1]="_" then the_new_letter:=parse(the_letter_as_string[2..]); RETURN(subs(the_letter=the_new_letter,the_sum)); else RETURN(the_sum); fi; end proc; if not(has(expr,Sum)) then RETURN(expr); else RETURN(subsindets( expr, 'specfunc( anything, Sum )', f->fix_it(f))); fi; end proc; sol:=dsolve(diff(y(x),x$2)+diff(y(x),x)+y(x)=0,y(x),'formal_series'): sol:=fix_summation_index(sol); y(x) = _C1*Sum((-1/2 - sqrt(3)*I/2)^n*x^n/n!, n = 0 .. infinity) + _C2*Sum((-1/2 + sqrt(3)*I/2)^n*x^n/n!, n = 0 .. infinity)
The latex now is
\[ y \left (x \right ) = \textit {\_C1} \left (\sum _{n =0}^{\infty }\frac {\left (-\frac {1}{2}+\frac {i \sqrt {3}}{2}\right )^{n} x^{n}}{n !}\right )+\textit {\_C2} \left (\sum _{n =0}^{\infty }\frac {\left (-\frac {1}{2}-\frac {i \sqrt {3}}{2}\right )^{n} x^{n}}{n !}\right ) \]