Posts Tagged ‘symbolic-math’

Convincing maple to solve an ODE with Neumann conditions at a symbolic valued location

Friday, November 17th, 2017

I can use maple to solve a 1D second-order ODE with Dirichlet boundary conditions at symbolic-valued locations:

# Z'' = 0, Z(a)=0, Z(b) = 1
dsolve({diff(Z(r),r,r) = 0,Z(a)=0,Z(b)=1}); 

This correctly returns

                  r       a
       Z(r) = - ----- + -----
                a - b   a - b

I can also easily convince maple to solve this ODE with some Neumann (normal derivative) boundary conditions at at fixed-value, numeric location:

# Z'' = 0, Z(a) = 1, Z'(0) = 0
dsolve({diff(Z(r),r,r) = 0,Z(a)=1,eval(diff(Z(r),r),r=0)=0});


                 Z(r) = 1

But if I try naively to use a Neumann condition at a symbolic value location

# Z'' = 0, Z(a) = 1, Z'(b) = 0
dsolve({diff(Z(r),r,r) = 0,Z(a)=1,eval(diff(Z(r),r),r=b)=0});

then I get an error:

Error, (in dsolve) found differentiated functions with same name but depending on different arguments in the given DE system: {Z(b), Z(r)}

After a long hunt, I found the solution. dsolve takes an optional second argument that can tell it what the dependent variable actually is. So the correct call is:

# Z'' = 0, Z(a) = 1, Z'(b) = 0
dsolve({diff(Z(r),r,r) = 0,Z(a)=1,eval(diff(Z(r),r),r=b)=0});

and this gives the correct answer

                 Z(r) = 1

Minimize a quadratic energy with linear equality constraints symbolically in maple

Wednesday, August 21st, 2013

I’m working with Maple to fit a polynomial to specific values and derivatives so that it minimized a quadratic energy. I tried the Optimization toolbox/package and came up with this:

p := (t) -> c3*t^3 + c2*t^2;
C1 := eval(p(t),t=1)-1;
C2 := eval(diff(p(t),t),t=1);
m := Minimize(int(diff(p(t),t)^2,t=0..1),{C1=0,C2=0});

but this seems to give me a numerically optimized solution:

m := [1.19999999999997, [c2 = 3., c3 = -2.]]

Rather I want a symbolically optimized solution. So far, I can do this by explicitly computing the Lagrangian and finding the saddle point:

Lambda := int(diff(p(t),t)^2,t=0..1) + lambda_1*(C1) + lambda_2*(C2);
s := solve(convert(Gradient(Lambda,[c3,c2,lambda_1,lambda_2]),list),[c3,c2,lambda_1,lambda_2]);

This gives me the exact solution:

[[c3 = -2, c2 = 3, lambda_1 = -12/5, lambda_2 = 1/5]]