Fwd: [maxima - Help] FEM Weak form & Maxima (Integration by Parts)



hello,

i'm forwarding this message which was posted in the
maxima forum at SF.
(https://sourceforge.net/forum/message.php?msg_id=3671925)
please direct followups to the mailing list instead of the forum.
to ck: you'll probably need to register in order to post
to the mailing list. sorry for the confusion.
see http://maxima.sf.net/maximalist.html .

robert dodier

--- "SourceForge.net" <noreply at sourceforge.net> wrote:

> To: noreply at sourceforge.net
> From: "SourceForge.net" <noreply at sourceforge.net>
> Subject: [maxima - Help] FEM Weak form & Maxima (Integration by
> Parts)
> Date: Wed, 05 Apr 2006 20:01:09 -0700
> 
> 
> Read and respond to this message at: 
> https://sourceforge.net/forum/message.php?msg_id=3671925
> By: ck_
> 
> 1) I'd like to be able to generate the weak form of an ODE/PDE
> entirely in Maxima.
> I'm a good part of the way there, but the partial integration by
> parts step
> is a road block.  Ideas?
> 
> sf: diff(u(x),x,2) - 4.0*u(x) = 0;
> sf1: expand(v(x)*sf);
> term_list(X) := block([], makelist(part(X,1,i), i,1,length(X)));
> terms: term_list(sf1);
> map( lambda([X], integrate(X,x,0,4.0)), terms );
> 
> /* Why cant Macsyma do integration by parts? */
> 
> wf1:integrate(diff(u(x),x)*diff(v(x),x),x,0,4) 
> 	+ integrate(4.0*v(x)*u(x),x,0,4) - v(x)*diff(u(x),x);
> dwf: wf1, 'diff(u(x),x,1)='sum(u(j)*diff(N(j,x),x,1),j,1,NN)$
> dwf: dwf, 'diff(v(x),x,1)='sum(v(i)*diff(N(i,x),x,1),i,1,NN)$
> dwf: dwf, u(x)='sum(u(j)*N(j,x),j,1,NN)$
> dwf: dwf, v(x)='sum(v(i)*N(i,x),i,1,NN)$
> 
> 2) Once you have the weak form, it would be nice to have Maxima
> generate the
> stiffness matrices as well... but the following chunk fails.  What is
> wrong?
> 
> /* Defining Shape Functions for a Quadratic Element */
> N1(x):= ((x-x2)*(x-x3)) / ((x1-x2)*(x1-x3))$
> N2(x):= ((x-x1)*(x-x3)) / ((x2-x1)*(x2-x3))$
> N3(x):= ((x-x1)*(x-x2)) / ((x3-x1)*(x3-x2))$
> dN1(x) := diff(N1(x),x)$
> dN2(x) := diff(N2(x),x)$
> dN3(x) := diff(N3(x),x)$
> Nij_quad: [N1(x),N2(x),N3(x)]$
> dNij_quad: map(lambda([X],diff(X,x,1)), Nij_quad)$
> 
> /* Defining Map for Quadratic Element Local Basis 
>  *	x2z: x=(L/2)*z + x2$
>  */
> x2z: x=((x3-x1)/2)*z + x2$
> 
> /* Specific Kije for this ODE */
> kij1e_quad(i,j) := 'integrate(dNij_quad[i]*dNij_quad[j],x,x1,x3)$
> kij1e_quad_z(i,j) := changevar(kij1e_quad(i,j),x2z,z,x)$
> kij2e_quad(i,j) := 'integrate(4.0*Nij_quad[i]*Nij_quad[j],x,x1,x3)$
> kij2e_quad_z(i,j) := changevar(kij2e_quad(i,j),x2z,z,x)$
> NN:3$
> kij1e_q:
> matrix([kij1e_quad_z(1,1),kij1e_quad_z(1,2),kij1e_quad_z(1,3)],
> 		[kij1e_quad_z(2,1),kij1e_quad_z(2,2),kij1e_quad_z(2,3)],
> 		[kij1e_quad_z(3,1),kij1e_quad_z(3,2),kij1e_quad_z(3,3)])$
> kij1e_q: ev(kij1e_q, nouns)$
> kij2e_q:
> matrix([kij2e_quad_z(1,1),kij2e_quad_z(1,2),kij2e_quad_z(1,3)],
> 		[kij2e_quad_z(2,1),kij2e_quad_z(2,2),kij2e_quad_z(2,3)],
> 		[kij2e_quad_z(3,1),kij2e_quad_z(3,2),kij2e_quad_z(3,3)])$
> kij2e_q: ev(kij2e_q, nouns)$
> 
> thanks all,
> 
> ck
> 
>
______________________________________________________________________
> You are receiving this email because you elected to monitor this
> forum.
> To stop monitoring this forum, login to SourceForge.net and visit: 
> https://sourceforge.net/forum/unmonitor.php?forum_id=15104
> 


__________________________________________________
Do You Yahoo!?
Tired of spam?  Yahoo! Mail has the best spam protection around 
http://mail.yahoo.com