<Text-field layout="Heading 1" style="Heading 1">Uebung 5</Text-field>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: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: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:
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: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="Heading 2" style="Heading 2">Aufgabe 1</Text-field><Text-field layout="Heading 3" style="Heading 3">a)</Text-field>restart;h:=(x)->1/sqrt(x^2-4);c:=3;d:=5;f := 10^(-5);Int(h(x),x=c..d);approx_t:=trapez(h,c,d,f);approx_m:=mittel(h,c,d,f);approx_s:=simpson(h,c,d,f);genau:=evalf(int(h(x),x=c..d));fehler_t:=abs(approx_t-genau);fehler_m:=abs(approx_m-genau);fehler_s:=abs(approx_s-genau);<Text-field layout="Heading 3" style="Heading 3">b)</Text-field>f:=(x)->x^2*exp(-x);a:=0;b:=1;e:=10^(-5);Int(f(x),x=a..b);approx_t:=trapez(f,a,b,e);approx_m:=mittel(f,a,b,e);approx_s:=simpson(f,a,b,e);genau:=evalf(int(f(x),x=a..b));fehler_t:=abs(approx_t-genau);fehler_m:=abs(approx_m-genau);fehler_s:=abs(approx_s-genau);<Text-field layout="Heading 2" style="Heading 2">Aufgabe 2</Text-field><Text-field layout="Heading 2" style="Heading 2">Aufgabe 3</Text-field><Text-field layout="Heading 3" style="Heading 3">a)</Text-field>restart;f:=(x)->1/sqrt(x^2-4);a:=3;b:=5;Int(f(x),x=a..b);res := [seq(romb(f,a,b,k,1),k=1..8)];genau:=evalf(int(f(x),x=a..b));seq(abs(op(i,res)-genau), i=1..8);res := [seq(romb(f,a,b,k,k),k=1..8)];genau:=evalf(int(f(x),x=a..b));seq(abs(op(i,res)-genau), i=1..8);<Text-field layout="Heading 3" style="Heading 3">b)</Text-field>restart;f:=(x)->x^2*exp(-x);a:=0;b:=1;Int(f(x),x=a..b);res := [seq(romb(f,a,b,k,1),k=1..8)];genau:=evalf(int(f(x),x=a..b));seq(abs(op(i,res)-genau), i=1..8);res := [seq(romb(f,a,b,k,k),k=1..8)];genau:=evalf(int(f(x),x=a..b));seq(abs(op(i,res)-genau), i=1..8);<Text-field layout="Heading 2" style="Heading 2">Aufgabe 4</Text-field><Text-field layout="Heading 3" style="Heading 3">a)</Text-field>restart;(a) int(ln(x+1)/(1-x)^(1/5), x = (0 .. 1)) = int(ln(1-x)/(x+1)^(1/5), x = (-1 .. 0)) = int((ln(1-x)-P[4](x))/(x+1)^(1/5), x = (-1 .. 0))+int(P[4](x)/(x+1)^(1/5), x = (-1 .. 0)) f:=(x)->ln(x+1)/(1-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) );a:=0;b:=1;n:=6;Int(f(x),x=a..b);genau := evalf( int(f(x),x=a..b) );approx := simpson_mit_n(g,6,-1,0) + evalf(int(p4(x)/(1+x)^(1/5), x=-1..0));abs(approx-genau);<Text-field layout="Heading 2" style="Heading 2">Aufgabe 5</Text-field>restart;with(plots, listplot):Digits := 10:g:=x->exp(-x^2/2)/sqrt(2*Pi);f:=x->simpson_mit_n(g,64,0,x) - c;df:=x->g(x);i:=1;e:=10^(-5);<Text-field layout="Heading 3" style="Heading 3">a)</Text-field>c:=0.481;genau := solve( int( g(t), t=0..x) = c, x );i:=1:
res:=[1]:
while abs(evalf(op(-1,res))-genau) > e do
printf("i = %d\n",i);
printf("fehler = %f\n", abs(evalf(op(-1,res))-genau));
tmp:=op(-1,res)-f(op(-1,res))/df(op(-1,res));
res:=[op(res),evalf(tmp)];
print(res);
i := i+1;
end do:plot(int(g(t),t=0..k), k=0..nops(res));<Text-field layout="Heading 3" style="Heading 3">b)</Text-field>c:=0.499;genau := solve( int( g(t), t=0..x) = c, x );i:=1;
res:=[1];
while abs(evalf(op(-1,res))-genau) > e do
printf("i = %d\n",i);
printf("fehler = %f\n", abs(evalf(op(-1,res))-genau));
tmp:=op(-1,res)-f(op(-1,res))/df(op(-1,res));
res:=[op(res),evalf(tmp)];
print(res);
i := i+1;
end do:plot(int(g(t),t=0..k), k=0..nops(res));<Text-field layout="Heading 3" style="Heading 3">c)</Text-field>c:=0.49999;genau := solve( int( g(t), t=0..x) = c, x );i:=1;
res:=[1];
while abs(evalf(op(-1,res))-genau) > e do
printf("i = %d\n",i);
printf("fehler = %f\n", abs(evalf(op(-1,res))-genau));
tmp:=op(-1,res)-f(op(-1,res))/df(op(-1,res));
res:=[op(res),evalf(tmp)];
print(res);
i := i+1;
end do:plot(int(g(t),t=0..k), k=0..nops(res));