<Text-field layout="Heading 1" style="Heading 1">Uebung 7</Text-field><Text-field layout="Heading 2" style="Heading 2">Aufgabe 1</Text-field>restart;with(plots,listplot):with(plots,display):#t:=(i) -> a + i*h;runge4 := proc(f1, f2, a, b, h, alpha1, alpha2)
local i, n, au1, au2, k11, k12, k21, k22, k31, k32, k41, k42, tmp1, tmp2;
n := (b-a)/h;
au1 := [alpha1];
au2 := [alpha2];
for i from 0 to n-1 do
k11 := evalf( h*f1(a+i*h, op(-1,au1), op(-1,au2)));
k12 := evalf( h*f2(a+i*h, op(-1,au1), op(-1,au2)));
k21 := evalf( h*f1( a+(i+1/2)*h,op(-1,au1)+k11/2, op(-1,au2)+k12/2));
k22 := evalf( h*f2( a+(i+1/2)*h,op(-1,au1)+k11/2, op(-1,au2)+k12/2));
k31 := evalf( h*f1( a+(i+1/2)*h,op(-1,au1)+k21/2, op(-1,au2)+k22/2));
k32 := evalf( h*f2( a+(i+1/2)*h,op(-1,au1)+k21/2, op(-1,au2)+k22/2));
k41 := evalf( h*f1( a+(i+1)*h,op(-1,au1)+k31, op(-1,au2)+k32));
k42 := evalf( h*f2( a+(i+1)*h,op(-1,au1)+k31, op(-1,au2)+k32));
tmp1:=op(-1,au1)+(1/6)*(k11+2*k21+2*k31+k41);
tmp2:=op(-1,au2)+(1/6)*(k12+2*k22+2*k32+k42);
au1:=[op(au1),tmp1];
au2:=[op(au2),tmp2];
end do;
return [au1,au2];
end proc:c1:=0.05;
c2:=7.5;################################################################
# y' = u2;
# y = u1;
#
# Gleichung umformen nach u3 und ersetzten
################################################################f1 := (t,u1,u2) -> u2;
f2 := (t,u1,u2) -> ;a:=0;b:=100;h:=(b-a)/3000;n:=(b-a)/h;au := runge4(f1, f2,a,b,h,1,0);plist := []:
for i from 1 to n do
tmp := [op(i,au[1]), op(i,au[2])];
plist := [op(plist), tmp];
end do:listplot(plist);<Text-field layout="Heading 2" style="Heading 2">Aufgabe 2</Text-field>restart;with(plots,listplot):with(plots,display):#t:=(i) -> a + i*h;runge4 := proc(f1, f2, a, b, h, alpha1, alpha2)
local i, n, au1, au2, k11, k12, k21, k22, k31, k32, k41, k42, tmp1, tmp2;
n := (b-a)/h;
au1 := [alpha1];
au2 := [alpha2];
for i from 0 to n-1 do
k11 := evalf( h*f1(a+i*h, op(-1,au1), op(-1,au2)));
k12 := evalf( h*f2(a+i*h, op(-1,au1), op(-1,au2)));
k21 := evalf( h*f1( a+(i+1/2)*h,op(-1,au1)+k11/2, op(-1,au2)+k12/2));
k22 := evalf( h*f2( a+(i+1/2)*h,op(-1,au1)+k11/2, op(-1,au2)+k12/2));
k31 := evalf( h*f1( a+(i+1/2)*h,op(-1,au1)+k21/2, op(-1,au2)+k22/2));
k32 := evalf( h*f2( a+(i+1/2)*h,op(-1,au1)+k21/2, op(-1,au2)+k22/2));
k41 := evalf( h*f1( a+(i+1)*h,op(-1,au1)+k31, op(-1,au2)+k32));
k42 := evalf( h*f2( a+(i+1)*h,op(-1,au1)+k31, op(-1,au2)+k32));
tmp1:=op(-1,au1)+(1/6)*(k11+2*k21+2*k31+k41);
tmp2:=op(-1,au2)+(1/6)*(k12+2*k22+2*k32+k42);
au1:=[op(au1),tmp1];
au2:=[op(au2),tmp2];
end do;
return [au1,au2];
end proc:c1 := 0.0008: c2 := 0.4:
c3 := 3.0: c4 := 0.002:
f1 := (t,x,y) -> x*(c1*y-c2);
f2 := (t,x,y) -> y*(c3-c4*x);a:=0;b:=100;h:=(b-a)/1000;n:=(b-a)/h;au := runge4(f1, f2,a,b,h,1000,2000);plist := []:
for i from 1 to n do
tmp := [op(i,au[1]), op(i,au[2])];
plist := [op(plist), tmp];
end do:listplot(plist);# stabile Loesungsolve( {x*(c1*y-c2) = 0, y*(c3-c4*x) = 0}, {x,y} );au := runge4(f1, f2,a,b,h,1500, 500);plist := []:
for i from 1 to n do
tmp := [op(i,au[1]), op(i,au[2])];
plist := [op(plist), tmp];
end do:listplot(plist);<Text-field layout="Heading 2" style="Heading 2">Aufgabe 3</Text-field>restart;with(plots,listplot);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:f := x -> sqrt(1+cos(x)^2):g := x -> simpson_mit_n(f,32,0,x):plot(g(x),x=1..50);# Lineare Interpolation zwischen x1=31 und x2=32 f\374r x=31.5
#y1 := g(31);
y2 := g(32);
lp := (y2-y1)*0.5+y1;
g(31.5);(g(32)+g(31))/2;# Vergleich Maple-intern
#l := [seq( evalf( int( f(x),x=0..k)), k=1..50)]:l;listplot(l);op(31,l);op(32,l);