\334bung 05 - 03.05.04<Text-field layout="_pstyle2" style="_cstyle2">Funktionen</Text-field><Text-field layout="_pstyle3" style="_cstyle3">List Tools</Text-field>Length := proc(ll) return nops(ll); end proc:
Last := proc(ll) return op(-1, ll); end proc:
Append := proc(ll, x) return [op(ll), x]; end proc:
Prepend := proc(ll, x) return [x, op(ll)]; end proc:
Member := proc(ll, p) return op(p, ll); end proc:
Take := proc(ll, n)
if n > 0 then return [ op(1..n, ll) ]
elif n < 0 then return [ op(n..-1, ll) ]
end if;
end proc:
Total := proc(ll)
local i;
return sum(Member(ll,i),i=1..Length(ll));
end proc:
Count := proc(ll, x)
local i,r;
r := 0;
for i from 1 by 1 to Length(ll) do
if Member(ll,i) = x then r := r+1 end if;
end do;
return r;
end proc:
Times := proc(l1, l2)
local i, r;
r := [];
if Length(l1) = Length(l2) then
for i from 1 to Length(l1) do
r := Append(r, Member(l1, i)*Member(l2, i));
end do:
end if;
return r;
end proc:with(plots, listplot):<Text-field layout="_pstyle3" style="_cstyle3">Einfache Quadraturen (Trapez, Mittelpunkt, Simpson)</Text-field># trapez: integriert fnc im Intervall [a,b] mit Schrittweite n nach der
# Trapezformel
##trapez := proc(fnc,n,a,b)
#end proc:trapez_n := proc(f,a,b,e)
local df2, res;
df2:=eval(diff(f(x),x$2),x=a);
res:=evalf(sqrt((b-a)^a*df2/(12*e)));
return round(res);
end proc:
trapez := proc(f,a,b,e)
local h, res, df2, n;
n:=trapez_n(f,a,b,e);
h:=(b-a)/n;
df2:=eval(diff(f(x),x$2),x=a);
#res:=h/2*(f(a)+f(b)+2*sum(f(a+j*h),j=1..n-1))-((b-a)*h^2)/12*df2;
res:=h/2*(f(a)+f(b)+2*sum(f(a+j*h),j=1..n-1));
return evalf(res);
end proc:
trapez_mit_n := proc(f,n,a,b)
local h, res, df2;
#n:=trapez_n(f,a,b,e);
h:=(b-a)/n;
df2:=eval(diff(f(x),x$2),x=a);
#res:=h/2*(f(a)+f(b)+2*sum(f(a+j*h),j=1..n-1))-((b-a)*h^2)/12*df2;
res:=h/2*(f(a)+f(b)+2*sum(f(a+j*h),j=1..n-1));
return evalf(res);
end proc:
# mitte: integriert fnc im Intervall [a,b] mit Schrittweite 2n nach der
# Mittelpunktsformel
#mitte := proc(fnc,n,a,b)
end proc:
mittel_n := proc(f,a,b,e)
local df2, res;
df2:=eval(diff(f(x),x$2),x=a);
res:=evalf(sqrt((b-a)^a*df2/(24*e)));
return round(res);
end proc:
mittel := proc(f,a,b,e)
local h, res, df2, n, sum, j;
n:=mittel_n(f,a,b,e);
h:=(b-a)/(n*2);
df2:=eval(diff(f(x),x$2),x=a);
sum:=0;
for j from 1 by 2 to 2*n-1 do
sum:=sum+f(a+j*h);
end do;
res := 2*h*sum;
return evalf(res);
end proc:
mittel_mit_n := proc(f,n,a,b)
local h, res, df2, sum, j;
h:=(b-a)/(n*2);
df2:=eval(diff(f(x),x$2),x=a);
sum:=0;
for j from 1 by 2 to 2*n-1 do
sum:=sum+f(a+j*h);
end do;
res := 2*h*sum;
return evalf(res);
end proc:# simpson: integriert fnc im Intervall [a,b] mit Schrittweite 2n nach der
# Simpson'schen Regel
#simpson := proc(fnc,n,a,b)
end proc:simpson_n := proc(f,a,b,e)
local df4, res;
df4 := eval(diff(f(x),x$4),x=a);
res:=evalf(((b-a)^b*df4/(16*180*e))^(1/4));
return round(res);
end proc:
simpson := proc(f,a,b,e)
local h, res, df4, n, sum1, sum2, j;
n:=simpson_n(f,a,b,e);
h:=(b-a)/(n*2);
df4 := eval(diff(f(x),x$4),x=a);
sum1:=0;
for j from 1 by 2 to 2*n-1 do
sum1:=sum1+f(a+j*h);
end do;
sum2:=0;
for j from 2 by 2 to 2*n-2 do
sum2:=sum2+f(a+j*h);
end do;
res:=h/3*(f(a)+f(b)+4*sum1+2*sum2);
#res:=h/3*(f(a)+f(b)+4*sum(f(a+j*h),j=1..2*n-1)+2*sum(f(a+j*h),j=2..2*n-2))-((b-a)*h^4)/180*df4;
return evalf(res);
end proc:
simpson_mit_n:=proc(f,n,a,b)
local h, res, df4, sum1, sum2, j;
h:=(b-a)/(n*2);
df4 := eval(diff(f(x),x$4),x=a);
sum1:=0;
for j from 1 by 2 to 2*n-1 do
sum1:=sum1+f(a+j*h);
end do;
sum2:=0;
for j from 2 by 2 to 2*n-2 do
sum2:=sum2+f(a+j*h);
end do;
res:=h/3*(f(a)+f(b)+4*sum1+2*sum2);
#res:=h/3*(f(a)+f(b)+4*sum(f(a+j*h),j=1..2*n-1)+2*sum(f(a+j*h),j=2..2*n-2))-((b-a)*h^4)/180*df4;
return evalf(res);
end proc:<Text-field layout="_pstyle3" style="_cstyle3">Gauss'sche Quadratur</Text-field>gauss := proc(f,n,a,b)
local r, c, res;
r:=[[0.5773502692, -0.5773502692], [0.7745966692, 0.0000000000, -0.7745966692], [0.8611363116, 0.3399810436, -0.3399810436, -0.8611363116], [0.9061798459, 0.5384693101, 0.00000000000, -0.5384693101, -0.9061798459]];
c:=[[1.0000000000, 1.0000000000], [0.5555555556, 0.8888888889, 0.5555555556], [0.3478548451, 0.6521451549, 0.6521451549, 0.3478548451], [0.2369268850, 0.4786286705, 0.5688888889, 0.4786286705, 0.2369268850]];
res:=(b-a)/2*sum(c[n-1][i]*f(((b-a)*r[n-1][i]+a+b)/2),i=1..n);
return res;
end proc:<Text-field layout="_pstyle3" style="_cstyle3">Romberg-Integration</Text-field># romb: Integriert nach Romberg-Schema
# romb := proc(f,a,b,k,j)
local res, h;
h:=(j)->(b-a)/(2^(j-1));
if(j = 1) then
if(k = 1) then
res:=(f(a)+f(b))*(b-a)/2;
else
res:=1/2*(romb(f,a,b,k-1,j)+h(k-1)*sum(f(a+(2*i-1)*h(k)),i=1..2^(k-2)));
end if;
else
res:=romb(f,a,b,k,j-1)+(romb(f,a,b,k,j-1)-romb(f,a,b,k-1,j-1))/(4^(j-1)-1);
end if;
return evalf(res);
end proc:<Text-field layout="_pstyle2" style="_cstyle2">Aufgabe 1</Text-field># Anzahl Schritte abh\344ngig von der Genauigkeit
#nm := (a,b,df2,e) -> evalf( sqrt( (b-a)^3*df2/(24*e) ) ):
nt := (a,b,df2,e) -> evalf( sqrt( (b-a)^3*df2/(12*e) ) ):
ns := (a,b,df2,e) -> evalf( ( (b-a)^5*df4/(16*180*e) )^(1/4) ):# int(1/sqrt(x^2-4), x=3..5)
#g := x -> 1/sqrt(x^2-4):
a := 3:
b := 5:
e := 10^(-5):tv := evalf( int(g(x), x=a..b) );plot(diff(g(x),x$2),x=a..b);df2 := eval(diff(g(x),x$2),x=3):
df4 := eval(diff(g(x),x$4),x=3):[ nm(a,b,df2,e), nt(a,b,df2,e), ns(a,b,df4,e) ];#[ mitte(g,115,a,b), trapez(g,162,a,b), simpson(g,8,a,b) ];
simpson_mit_n(g,8,a,b):
simpson1:=evalf(%);
trapez_mit_n(g,162,a,b):
trapez1:=evalf(%);
mittel1:=evalf(mittel(g,a,b,e));g := x -> sqrt(1+x^2):
a := 1:
b := Pi:
e := 10^(-4):
exakt := evalf( int(g(x), x=a..b) );
df2 := eval(diff(g(x),x$2),x=3):
df4 := eval(diff(g(x),x$4),x=3):
#nss:= (a,b,df2,e) -> evalf( ( (b-a)^5*df4/(16*180*e) )^(1/4) )
evalf(ns(a,Pi,df4,e));
smpsn:=simpson(g,a,b,e);# int(x^2*exp(-x), x=0..1)
#g := x -> x^2*exp(-x):
a := 0:
b := 1:
e := 10^(-5):tv := evalf( int(g(x), x=a..b) );plot(diff(g(x),x$2),x=a..b);df2 := eval(diff(g(x),x$2),x=0):
df4 := eval(diff(g(x),x$4),x=0):[ nm(a,b,df2,e), nt(a,b,df2,e), ns(a,b,df4,e) ];[ mitte(g,92,a,b), trapez(g,130,a,b), simpson(g,5,a,b) ];<Text-field layout="_pstyle2" style="_cstyle2">Aufgabe 2</Text-field># int(1/sqrt(x^2-4), x=3..5)
# g := x -> 1/sqrt(x^2-4):
a := 3:
b := 5:tv := evalf( int(g(x),x=a..b) ):tv;
res := [seq(gauss(g,k,a,b),k=2..5)];seq(abs( Member(res,i-1)-tv), i=2..5);# int(x^2*exp(-x), x=0..1)
#g := x -> x^2*exp(-x):
a := 0:
b := 1:tv := evalf( int(g(x),x=a..b) ):tv;
res := [seq(gauss(g,k,a,b),k=2..5)];seq(abs( Member(res,i-1)-tv), i=2..5);<Text-field layout="_pstyle2" style="_cstyle2">Aufgabe 3</Text-field># int(1/sqrt(x^2-4), x=3..5)
# g := x -> 1/sqrt(x^2-4):
a := 3:
b := 5:tv := evalf( int(g(x),x=a..b) ):# Erste Ordnung
#tv;
res := [seq(romb(g,a,b,k,1),k=1..8)];err := [ seq(abs( Member(res,i)-tv), i=1..8) ];# Beschleunigt
#res := [seq(romb(g,a,b,k,k),k=1..8)];err := [ seq(abs( Member(res,i)-tv), i=1..8) ];# int(x^2*exp(-x), x=0..1)
# g := x -> x^2*exp(-x):
a := 0:
b := 1:tv := evalf( int(g(x),x=a..b) ):# Erste Ordnung
#tv;
res := [seq(romb(g,a,b,k,1),k=1..8)];err := [ seq(abs( Member(res,i)-tv), i=1..8) ];# Beschleunigt
#res := [seq(romb(g,a,b,k,k),k=1..8)];err := [ seq(abs( Member(res,i)-tv), i=1..8) ];<Text-field layout="_pstyle2" style="_cstyle2">Aufgabe 4</Text-field>(a) NiMvLSUkaW50RzYkKiYtJSNsbkc2IywmJSJ4RyIiIkYtRi1GLSksJkYtRi1GLCEiIiomRi1GLSIiJkYwRjAvRiw7IiIhRi0tRiU2JComLUYpNiNGL0YtKUYrRjFGMC9GLDssJEYtRjBGNQ== = NiMsJi0lJGludEc2JComLCYtJSNsbkc2IywmIiIiRi0lInhHISIiRi0tJiUiUEc2IyIiJTYjRi5GL0YtKSwmRi5GLUYtRi0qJkYtRi0iIiZGL0YvL0YuOywkRi1GLyIiIUYtLUYlNiQqJkYwRi1GNkYvRjpGLQ== # Genauer Wert
#tv := evalf( int(ln(x+1)/(1-x)^(1/5),x=0..1) );evalf( int(ln(1-x)/(x+1)^(1/5),x=-1..0) );# Approximation
#taylor(ln(1-x), x=-1, 5);p4 := x -> ln(2)-1/2*(x+1)-1/8*(x+1)^2-1/24*(x+1)^3-1/64*(x+1)^4:g := x -> piecewise( x=-1, 0, x<>-1, (ln(1-x)-p4(x))/(1+x)^(1/5) ):approx := simpson_mit_n(g,6,-1,0)+evalf( int(p4(x)/(1+x)^(1/5), x=-1..0) );abs(approx-tv);# Aufgabe 2)exakt := evalf( int(exp(2*x)/x^(2/5),x=0..1) );f:=x->exp(2*x);p4 := x -> ln(2)-1/2*(x+1)-1/8*(x+1)^2-1/24*(x+1)^3-1/64*(x+1)^4;taylor(f(x), x=0, 5);p4:=x->1+2*x^2+2*x^2+4/3*x^3+2/3*x^4;g := x -> piecewise( x=0, 1, x<>0, (f(x)-p4(x))/x^(2/5) );
evalf(%);approx := simpson_mit_n(g,6,0,1)+evalf( int(p4(x)/x^(2/5), x=0..1) );(b) NiMvLSUkaW50RzYkKiYtJSRjb3NHNiMqJiIiIkYsJSJ4RyEiIkYsKiRGLSIiJEYuL0YtO0YsJSlpbmZpbml0eUctRiU2JComRi1GLC1GKTYjRi1GLC9GLTsiIiFGLA== ( In der Aufgabenstellung heisst es irrt\374mlicherweise NiMtJSRpbnRHNiQqJi0lJGNvc0c2IyUieEciIiIqJEYqIiIkISIiL0YqO0YrJSlpbmZpbml0eUc= ) # Genauer Wert
#tv := evalf( int(cos(1/x)/x^3, x=1..infinity) );evalf( int(x*cos(x),x=0..1) );# Approximation
#g := x -> x*cos(x):approx := simpson(g,6,0,1);abs(approx-tv);