Method of Weighted Residuals Commas Lecture "Discretization Methods", held by Dr.-Ing. A. Schmidt Exercise 3 Maple file created by Stefan Spreng in Dec 2009

restart:with(Optimization): numerical values that will be used later to plot the results (not used for the computations): lvalue:=1: # length qvalue:=1: # magnitude of the load EIvalue:=1: # bending stiffness define a function exactsolution(x) for the exact solution as a reference: exactsolution:=proc(x) q*l^4/(360*EI) * (3*(x/l)^5 - 10*(x/l)^3 + 7*(x/l)): end proc: exactsolution(x); 1. Write down the DE and the BC differential equation: defined in the form DE=0 (with u4 being the 4th derivative of u) DE:=EI*u4-q/l*x; specify number of BC and define BC: n_BC:=4; BC:=[[x=0, u=0], [x=l, u=0], [x=0, u2=0], [x=l, u2=0]]; 2. Choose number n of free parameters c_i define the number n of free parameters: n:=4; define the list of free parameters c c:=[seq(c||i, i=1..n)]; 3. Choose n appropriate shape functions define n shape functions s: s:=[seq(sin(i*Pi*x/l), i=1..n)]: for i from 1 to n do print(s||i= s[i]); od; plot([seq(eval(s[i], l=lvalue), i=1..n)], x=0..lvalue, legend=[seq(cat("s",j), j=1..n)]); 4. Assemble the trial function u u = sum (c_i * s_i ): u:=sum(c[j]*s[j],j=1..n); find the necessary derivatives of the trial function: u1:=combine(diff(u,x),trig); u2:=combine(diff(u1,x),trig); u3:=combine(diff(u2,x),trig); u4:=combine(diff(u3,x),trig); we use an interior method: check whether the boundary conditions are fulfilled for i from 1 to n_BC do if (eval(subs(BC[i][1], BC[i][2]))) then print("BC "||i||" ok"); else print("BC "||i||" NOT ok"); end if; od; 5. Define the Residual define the residual as the differential equation in terms of the trial function, i.e. as a function of the free parameters: R:=DE; now we need a way to determine the free parameters c_1 to c_n. 4 different methods to do this will be shown:

The method of Subdomains is a method to find values for the free parameters in the trial function unassign('SDcollpoints', 'SDeqns', 'SDc', 'SDu', 'SDfuncu', 'SDfuncuf', 'approx', 'exact', 'l', 'q', 'EI', 'x1'); 6a. Choose weight functions and create n equations for the free parameters define n subdomains of equal length: SDsubdom:=[seq(l/n*i, i=0..n)]; create n equations for the free parameters: for i from 1 to n do SDeqns[i]:=expand(int(R,x=SDsubdom[i]..SDsubdom[i+1]))=0; od; 7a. Solve the system of equations for the free parameters SDc:=(solve({seq(SDeqns[i],i=1..n)},c))[1]; # solve returns a list of lists -> use only the first entry 8a. Rewrite the trial function with the solution for the free parameters This is the approximate solution to the DE found with this method SDu:=subs(SDc,u); 9a. Check the Results define functions SDfuncu(x) and SDfuncuf(x) that return the result obtained from the trial function at point x as an exact value or a numerical value respectively SDfuncu:=proc(xx) sort(simplify(eval(subs(x=xx,SDu)))); end proc: SDfuncuf:=proc(xx) evalf(subs(x=xx,SDu)); end proc: print the result at the point x1: x1:=l/2; # should be given in terms of l approx:=SDfuncu(x1)=SDfuncuf(x1); exact:=exactsolution(x1)=evalf(exactsolution(x1)); plot the results after assuming numerical values for l, q, and EI: l:=lvalue: q:=qvalue: EI:=EIvalue: plot([-SDfuncu(x),-exactsolution(x)], x=0..l, title = "comparison of exact and approximated results", labels=["length", "deflection"], gridlines=true, tickmarks=[10, default], legend=["approx.","exact"]); SDs:=seq(subs(SDc[i], -c[i]*s[i]), i=1..n): SDlegend:=seq(cat("j = ",i), i=1..n): plot([-SDfuncu(x), SDs], x=0..l, title = "how do the n different shape functions add up?", labels=["length", "c_j * u_j"], gridlines=true, tickmarks=[10, default], legend=["sum",SDlegend]); Error analysis: absolute error SDerr:=proc(x) evalf((exactsolution(x)-SDfuncu(x))); end proc: plot(SDerr(x), x=0..l, labels=["length", "absolute error"], gridlines=true); Maximum value of absolute error SDexterr:=Maximize(abs(SDerr(x)),x=0..l); Relative error at this location SDextrelerr:=SDexterr[1]/eval(exactsolution(x),SDexterr[2]); unassign('approx', 'exact', 'l', 'q', 'EI', 'x1'): sort([anames(user)]):

The Collocation Method is a method to find values for the free parameters in the trial function unassign('CMcollpoints', 'CMeqns', 'CMc', 'CMu', 'CMfuncu', 'CMfuncuf', 'approx', 'exact', 'l', 'q', 'EI', 'x1'); 6b. Choose weight functions and create n equations for the free parameters define n collocation points in the domain l: (you can do that either automatically or by hand. If done automatically, the collocation points have equal distances. If done manually, you have to take care that you adjust the number of collocation points to the numberr n of free parameters!) CMcollpoints:=[seq(l/(n+1)*i, i=1..n)]; # this line creates n collocation points with equal distance #CMcollpoints:=[1/2*l, 3/4*l]; # here, the collocation points can be entered manually create n equations for the free parameters: for i from 1 to n do CMeqns[i]:=expand(subs(x=CMcollpoints[i], R))=0; od; 7b. Solve the system of equations for the free parameters CMc:=(solve({seq(CMeqns[i],i=1..n)},c))[1]; # solve returns a list of lists -> use only the first entry 8b. Rewrite the trial function with the solution for the free parameters This is the approximate solution to the DE found with this method CMu:=subs(CMc,u); 9b. Check the Results define functions CMfuncu(x) and CMfuncuf(x) that return the result obtained from the trial function at point x as an exact value or a numerical value respectively CMfuncu:=proc(xx) sort(simplify(eval(subs(x=xx,CMu)))); end proc: CMfuncuf:=proc(xx) evalf(subs(x=xx,CMu)); end proc: print the result at the point x1: x1:=l/2; # should be given in terms of l approx:=CMfuncu(x1)=CMfuncuf(x1); exact:=exactsolution(x1)=evalf(exactsolution(x1)); plot the results after assuming numerical values for l, q, and EI: l:=lvalue: q:=qvalue: EI:=EIvalue: plot([-CMfuncu(x), -exactsolution(x)], x=0..l, title = "comparison of exact and approximated results", labels=["length", "deflection"], gridlines=true, tickmarks=[10, default], legend=["approx.","exact"]); CMs:=seq(subs(CMc[i], -c[i]*s[i]), i=1..n): CMlegend:=seq(cat("j = ",i), i=1..n): plot([-CMfuncu(x), CMs], x=0..l, title = "how do the n different shape functions add up?", labels=["length", "c_j * s_j"], gridlines=true, tickmarks=[10, default], legend=["sum",CMlegend]); Error analysis: absolute error CMerr:=proc(x) evalf((exactsolution(x)-CMfuncu(x))); end proc: plot(CMerr(x), x=0..l, labels=["length", "absolute error"], gridlines=true); Maximum value of absolute error CMexterr:=Maximize(abs(CMerr(x)),x=0..l); Relative error at this location CMextrelerr:=CMexterr[1]/eval(exactsolution(x),CMexterr[2]); unassign('approx', 'exact', 'l', 'q', 'EI', 'x1'): sort([anames(user)]):

The Least Squares Method is a method to find values for the free parameters in the trial function unassign('LSweightfunctions', 'LSeqns', 'LSc', 'LSu', 'LSfuncu', 'LSfuncuf', 'approx', 'exact', 'l', 'q', 'EI', 'x1'); 6c. Choose weight functions and create n equations for the free parameters define n weight functions of the form dR / dc_j for i from 1 to n do LSweightfunctions[i]:=diff(R, c[i]); od; create n equations for the free parameters: for i from 1 to n do LSeqns[i]:=expand(int(LSweightfunctions[i]*R, x=0..l))=0; od; 7c. solve the system of equations for the free parameters LSc:=(solve({seq(LSeqns[i],i=1..n)},c))[1]; # solve returns a list of lists -> use only the first entry 8c. Rewrite the trial function with the solution for the free parameters This is the approximate solution to the DE found with this method LSu:=subs(LSc,u); 9c. Check the Results define functions LSfuncu(x) and LSfuncuf(x) that return the result obtained from the trial function at point x as an exact value or a numerical value respectively LSfuncu:=proc(xx) sort(simplify(eval(subs(x=xx,LSu)))); end proc: LSfuncuf:=proc(xx) evalf(subs(x=xx,LSu)); end proc: print the result at the point x1: x1:=l/2; # should be given in terms of l approx:=LSfuncu(x1)=LSfuncuf(x1); exact:=exactsolution(x1)=evalf(exactsolution(x1)); plot the results after assuming numerical values for l, q, and EI: l:=lvalue: q:=qvalue: EI:=EIvalue: plot([-LSfuncu(x),-exactsolution(x)], x=0..l, title = "comparison of exact and approximated results", labels=["length", "deflection"], gridlines=true, tickmarks=[10, default], legend=["approx.","exact"]); LSs:=seq(subs(LSc[i], -c[i]*s[i]), i=1..n): LSlegend:=seq(cat("j = ",i), i=1..n): plot([-LSfuncu(x), LSs], x=0..l, title = "how do the n different shape functions add up?", labels=["length", "c_j * u_j"], gridlines=true, tickmarks=[10, default], legend=["sum",LSlegend]); Error analysis: absolute error LSerr:=proc(x) evalf((exactsolution(x)-LSfuncu(x))); end proc: plot(LSerr(x), x=0..l, labels=["length", "absolute error"], gridlines=true); Maximum value of absolute error LSexterr:=Maximize(abs(LSerr(x)),x=0..l); Relative error at this location LSextrelerr:=LSexterr[1]/eval(exactsolution(x),LSexterr[2]); unassign('approx', 'exact', 'l', 'q', 'EI', 'x1'): sort([anames(user)]):

Galerkin's Method is a method to find values for the free parameters in the trial function unassign('GMweightfunctions', 'GMeqns', 'GMc', 'GMu', 'GMfuncu', 'GMfuncuf', 'approx', 'exact', 'l', 'q', 'EI', 'x1'); 6d. Choose weight functions and create n equations for the free parameters define n weight functions which are the same as the shape functions for i from 1 to n do GMweightfunctions[i]:=s[i]; od; create n equations for the free parameters: for i from 1 to n do GMeqns[i]:=(int(GMweightfunctions[i]*R, x=0..l))=0; od; 7d. Solve the system of equations for the free parameters GMc:=(solve({seq(GMeqns[i],i=1..n)},c))[1]; # solve returns a list of lists -> use only the first entry 8d. Rewrite the trial function with the solution for the free parameters This is the approximate solution to the DE found with this method GMu:=subs(GMc,u); 9d. Check the Results define functions GMfuncu(x) and GMfuncuf(x) that return the result obtained from the trial function at point x as an exact value or a numerical value respectively GMfuncu:=proc(xx) sort(simplify(eval(subs(x=xx,GMu)))); end proc: GMfuncuf:=proc(xx) evalf(subs(x=xx,GMu)); end proc: print the result at the point x1: x1:=l/2; # should be given in terms of l approx:=GMfuncu(x1)=GMfuncuf(x1); exact:=exactsolution(x1)=evalf(exactsolution(x1)); plot the results after assuming numerical values for l, q, and EI: l:=lvalue: q:=qvalue: EI:=EIvalue: plot([-GMfuncu(x),-exactsolution(x)], x=0..l, title = "comparison of exact and approximated results", labels=["length", "deflection"], gridlines=true, tickmarks=[10, default], legend=["approx.","exact"]); GMs:=seq(subs(GMc[i], -c[i]*s[i]), i=1..n): GMlegend:=seq(cat("j = ",i), i=1..n): plot([-GMfuncu(x), GMs], x=0..l, title = "how do the n different shape functions add up?", labels=["length", "c_j * u_j"], gridlines=true, tickmarks=[10, default], legend=["sum",GMlegend]); Error analysis: absolute error GMerr:=proc(x) evalf((exactsolution(x)-GMfuncu(x))); end proc: plot(GMerr(x), x=0..l, labels=["length", "absolute error"], gridlines=true); Maximum value of absolute error GMexterr:=Maximize(abs(GMerr(x)),x=0..l); Relative error at this location GMextrelerr:=GMexterr[1]/eval(exactsolution(x),GMexterr[2]); unassign('approx', 'exact', 'l', 'q', 'EI', 'x1'): sort([anames(user)]):

Graphical comparison of the different solutions: l:=lvalue: q:=qvalue: EI:=EIvalue: plot([-exactsolution(x),-SDfuncu(x),-CMfuncu(x),-LSfuncu(x),-GMfuncu(x)],x=0..l, legend=["exact", Subdomains, Collocation, LeastSquares, Galerkin]); Comparison of the absolute errors plot([SDerr(x),CMerr(x),LSerr(x),GMerr(x)],x=0..l, legend=[Subdomains, Collocation, LeastSquares, Galerkin]); Tabular overview of the different values for the free parameters: (Rigth click into the table an choose "Evaluate All" to obtain current values) NiNJK1N1YmRvbWFpbnNHNiI= NiNJLENvbGxvY2F0aW9uRzYi NiNJLUxlYXN0U3F1YXJlc0c2Ig== NiNJKUdhbGVya2luRzYi NiNJI2MxRzYi NiNJI2MyRzYi NiNJI2MzRzYi NiNJI2M0RzYi NiNJI2M1RzYi NiMvLCQqJiIlKzUiIiJJI2MxRzYiIiIiIiIiJCIrME8pRzEpISIq NiMvLCQqJiIlKzUiIiJJI2MyRzYiIiIiIiIiJCErRjZsPkQhIzU= Error, invalid subscript selector Error, invalid subscript selector Error, invalid subscript selector NiMvLCQqJiIlKzUiIiJJI2MxRzYiIiIiIiIiJCIrO3cxRmYhIio= NiMvLCQqJiIlKzUiIiJJI2MyRzYiIiIiIiIiJCErYGQhW0IiISM1 Error, invalid subscript selector Error, invalid subscript selector Error, invalid subscript selector NiMvLCQqJiIlKzUiIiJJI2MxRzYiIiIiIiIiJCIrI0dGYmAnISIq NiMvLCQqJiIlKzUiIiJJI2MyRzYiIiIiIiIiJCErd0FOVT8hIzU= Error, invalid subscript selector Error, invalid subscript selector Error, invalid subscript selector NiMvLCQqJiIlKzUiIiJJI2MxRzYiIiIiIiIiJCIrI0dGYmAnISIq NiMvLCQqJiIlKzUiIiJJI2MyRzYiIiIiIiIiJCErd0FOVT8hIzU= Error, invalid subscript selector Error, invalid subscript selector Error, invalid subscript selector NiMiJSs1 NiNJL3JlYWx0aXZlfmVycm9yRzYi NiMkIitHeU4nUSMhIzU= NiMkIitWRVxWJiohIzY= NiMkIitfJGVtOSIhIzY= NiMkIitfJGVtOSIhIzY=

JSFH