całka simpsona w pascalu

witam
mam problem z napisaniem programu w pascalu ktory bedzie liczyc calke metoda simpsona
calka wyglada tak:

mam taki programik ale daje zły wynik , poprawny jest : 0.513
a dokładność ma być równa 10 do potegi -4
///////////////////////////////////////////////////////////////////////
program sim;
uses crt;
{var
pom:real;}

function F(x:real):real;
begin
F:=x/1+(sqrt(ln(x)));
end;

function simp (m:integer; a,b:real):real; { nie wiem po co zmienna m }
var
i:integer;
x,rk,s:real;
begin
rk:=(b-a)/m;
s:=0.0;
for i:=1 to m do
begin
x:=a+rk/2+(i-1)rk;
s:=s+4f(x);
end;
for i:=1 to m-1 do
begin
x:=a+irk;
s:=s+2f(x)
end;
simp:=(rk/6)*(f(a)+s+f(b))
end;
{--------------------------------------------------------------------}
begin
clrscr;
writeln;
textcolor(15);
writeln;
writeln('Dana jest calka w przedziale od 2 do 3');
writeln;
writeln('Rozwiazanie za pomoca metody Simpsona: ');
writeln(simp(5,2.0,3.0)7);
readkey;
end.
/////////////////////////////////////////////////////////////////
byłbym wdzieczny jak by ktoś poprawił błędy i napisał jakis króciutki komentarz
zebym to zrozumiał i wiedział na przyszłość
pozdro

może jakieś info jak się taką całke liczy bo dawno to było :P

Proponuje napisac funkcje F w postaci :)

function F(x:real):real;
begin
  F:=1/(1+sqrt(ln(x)) );
end;

PS. Bawic sie latex

program sim;
uses crt;
{var
pom:real;}

function F(x:real):real;
begin
  F:=1/(1+sqrt(ln(x)));{to jest prawidowo zdefiniowana funkcja}
end;

function simp (m:integer; a,b:real):real; { m to ilosc podzialow przedzailu }
var
  i:integer;
  x,rk,s:real;
begin
  rk:=(b-a)/m; {tu jest potrzebne to m i powinno byc rzedu 100-200}
              { dzielisz przedzial (3-2) na m czesci}     
  s:=0.0;
  for i:=1 to m do
    begin
      x:=a+rk/2+(i-1)*rk;
      s:=s+4*f(x);
    end;
  for i:=1 to m-1 do
    begin
      x:=a+i*rk;
      s:=s+2*f(x)
    end;
  simp:=(rk/6)*(f(a)+s+f(b))
end;
{--------------------------------------------------------------------}
begin
  clrscr;
  writeln;
  textcolor(15);
  writeln;
  writeln('Dana jest calka w przedziale od 2 do 3');
  writeln;
  writeln('Rozwiazanie za pomoca metody Simpsona: ');
  writeln(simp(200,2.0,3.0):3:4);
  Readln; {program poczeka na wcisnięcie ENTER i bedziesz widzial wynik}
end.

siemka
StPan - wielkie dzieki za pomoc
juz teraz widze ze zle zdefiniowalem funkcje [sciana]
no juz teraz to prawie rozumiem ale bylbym wdzieczny jakbys jeszcze
mi dopisal chociaz dlaczego "m" powinno byc rzedu 100-200? bo jesli
wezme inna liczbe ( np 5 ; 500 ) to tez w sumie dziala , a tak z
matematycznego punktu to po co sa te podziały przedziału :| ?
jeszcze raz dzieki
pozdro

czesc, aby liczyc calke metoda sipmsona (czyli w zasadzie liczyc ja lopatologicznie bo z definicji czyli dzielisz pole pod krzywa na n prostokatow i sumujesz pola tych prostokatow) musisz ustalic na ile przedzialow podzielisz przedzial calkowania. Im wiecej tym lepsza dokladnosc ale dluzszy czas liczenia, no i ilosc ustala sie tez w zaleznosci do tego w jakim przedziale chcesz liczyc calke - jak jest to przedzial [-100,100] to wiadomo ze n=10 to za malo bo prostokaty beda mialy zbyt duza szerokosc i slabo beda przyblizaly faktyczne pole pod krzywa (mozna uzyc jeszcze metody trapezow jest to metoda b.podobna)

Pozdr. skalniak

tylko, ze metoda trapezow daje mniej dokladne wyniki.

bob999 napisał(a)

tylko, ze metoda trapezow daje mniej dokladne wyniki.

hmmm wydawalo mi sie ze jednak daje lepsze bo w polowie chociaz uwzgledniasz to pole ktorego nie uwzgledniasz przy uzyciu pol prostokatow

Pozdr. skalniak

Mylisz metody.

• prostokątów - rzędu 0
• trapezów - 1
• parabol, czyli Simpsona - 2
• wielomianami 3-go stopnia - 3
  ...

Są to dość marne metody.

Dobra jest np. metoda Romberga.

zapraszam do lektury:
http://www.i-lo.tarnow.pl/edu/inf/alg/calki/

macie jeszcze met. Romberg'a

function F(x:real):real;
begin
  F:=x/1+(sqrt(sin(x))); //definiujemy sobie dowolna funkcje

end;


PROCEDURE tform1.romb(lower, upper, tol: extended; VAR ans: extended);
{ numerical integration by the Romberg method }
{ function fx, name cannot be passed by Turbo Pascal}
VAR
	nx			: ARRAY[1..16] of Integer;
	t			: ARRAY[1..136] of Real;
	done,error		: Boolean;
	pieces,nt,i,ii,n,nn,
	l,ntra,k,m,j		: Integer ;
	delta_x,c,sum,fotom,x	: Real ;
BEGIN
  done := false;
  error := false;
  pieces := 1;
  nx[1] := 1;
  delta_x := (upper-lower)/pieces;
  c := (F(lower)+F(upper))*0.5;
  t[1] := delta_x*c;
  n := 1;
  nn := 2;
  sum := c;
  REPEAT
    n := n+1;
    fotom := 4.0;
    nx[n] := nn;
    pieces := pieces*2;
    l := pieces-1;
    delta_x := (upper-lower)/pieces;
    { compute trapezoidal sum for 2^(n-1)+1 points }
    FOR ii:=1 to (l+1) div 2 DO
      BEGIN
	i := ii*2-1;
	x := lower+i*delta_x;
	sum := sum+F(x)
      END;
    t[nn] := delta_x*sum;
    //memo1.Lines.Add(inttostr(pieces));
    ntra := nx[n-1];
    k := n-1;
    { compute n-th row of T array }
    FOR m:=1 to k DO
      BEGIN
	j := nn+m;
	nt := nx[n-1]+m-1;
	t[j] := (fotom*t[j-1]-t[nt])/(fotom-1.0);
	fotom := fotom*4.0
      END;
      //memo1.Lines.Add(inttostr(j));
    //writeln(j:4,t[j]);
    IF n>4 THEN
      BEGIN
	IF t[nn+1]<>0.0 THEN
	  IF (abs(t[ntra+1]-t[nn+1])<=abs(t[nn+1]*tol))
	    OR (abs(t[nn-1]-t[j])<=abs(t[j]*tol)) THEN
	      done := true
	  ELSE IF n>15 THEN
	    BEGIN
	      done := true;
	      error := true
	    END
	END;	{ if n>4 }
    nn := j+1
  UNTIL done;
  ans := t[j];  //ans zmienna przechowujaca wynik

END;

function fx, name cannot be passed by Turbo Pascal

Autor tego kodu był jakiś niedouczony.
Procedurę można przekazać jako parametr w dowolnym kompilatorze paskala.

kwadraturowiec napisał(a)

Mylisz metody.
Dobra jest np. metoda Romberga.

jedna z lepszych lecz dosc trudnych w porownaniu z wymienianymi jest metoda Gaussa-Legendrea :), nieziemsko szybka:)

Pozdr. skalniak

skalniak napisał(a)

jedna z lepszych lecz dosc trudnych w porownaniu z wymienianymi jest metoda Gaussa-Legendrea :), nieziemsko szybka:)

To jest metoda innego typu - węzły nie są równomiernie rozmieszczone,
lecz są to pierwiastki wielomianów Legendrea - i to jest pewien problem.
oto dwa takie wielomiany:
35x^{4 - 30x}2 + 3 = 0, otrzymamy cztery węzły
231x^{6 - 315x}3 + 105x^2 - 5 = 0, a tu 6

Można stablicować te pierwiastki oraz współczynniki i wyjdzie metoda,
chyba nawet, prostsza od Romberga.

Nie jest pewne czy metoda ta będzie szybsza.
Jeśli błąd jest zbyt duży, wszystkie obliczenia robimy od nowa - z większą liczbą przedziałów.
A w Romberga korzystamy z wcześniejszych obliczeń! :)

kwadraturo napisał(a)

skalniak napisał(a)

jedna z lepszych lecz dosc trudnych w porownaniu z wymienianymi jest metoda Gaussa-Legendrea :), nieziemsko szybka:)

To jest metoda innego typu - węzły nie są równomiernie rozmieszczone,
lecz są to pierwiastki wielomianów Legendrea - i to jest pewien problem.
oto dwa takie wielomiany:
35x^{4 - 30x}2 + 3 = 0, otrzymamy cztery węzły
231x^{6 - 315x}3 + 105x^2 - 5 = 0, a tu 6

Można stablicować te pierwiastki oraz współczynniki i wyjdzie metoda,
chyba nawet, prostsza od Romberga.

Nie jest pewne czy metoda ta będzie szybsza.
Jeśli błąd jest zbyt duży, wszystkie obliczenia robimy od nowa - z większą liczbą przedziałów.
A w Romberga korzystamy z wcześniejszych obliczeń! :)

nie musisz mi tlumaczyc znam metode ktora proponuje:) nie jest to problemem ze trzeba znajdowac pierwiastki wielomianu, mozna go zapisac w postaci szeregow (bo normalna postac jest rozniczkowa i nieprzyjazna do zaprogramowania) i obliczac dla dowolnego n gdzie n jest iloscia punktow - nalezy zaimplementowac rozwiazanie ukladu n-rownan z n -niewiadomymi, mozna nawet wszystko wyliczyc wczesniej np dla 1000 punktow, trzymac w pliku - dla duzej dokladnosci znacznie przyspieszy to obliczenia. :) Tej metody Romberga nie znam w ogóle.

Pozdr. skalniak

Ogólnie wiadomo, że kwadratury wysokiego rzędu są złym rozwiązaniem (zjawisko Rungego),
a kwadratura gaussa z tysiącem węzłów, to już chyba rekord świata.

Przedzał całkowania dzielimy na mniejsze, i w każdym stosujemy kwadraturę niskiego rzędu.

umiałby ktoś napisać program w Pascalu na obliczenie całki oznaczonej z (xlnx) w granicach 3 do 5 załóżmy metoda trapezów?

a) odkopujesz temat sprzed prawie 5 lat :D,
b) bym się strasznie zdziwił, gdyby to nie było dostępne na necie,
c) jeżeli ma działać to na formie "napisać na zaliczenie" to zapytanie do działu Praca.

A po co mam zakładać nowy temat? Akurat w necie nic takiego nie ma, poza tym potrzebuje jeszcze program na wyliczenie reakcji w belce, a tego to juz na pewno nie ma, ale nie za bardzo rozumiem, co masz na mysli w podpunkcie c)?

ok czyli rozumiem, że nikt mi nie pomoże?

Ja z matmy to jestem "noga", ale tak jak Tobie napisano poprzednio - nie masz wcale samodzielnego kodu, to
jak potrzebujesz mieć program na zaliczenie to napisz do działu Offtopic - Praca proponując tam kwotę za kod.
A to forum jest po to żeby z problemem pomóc kogoś naprowadzić czy poprawić, a nie odwalić robotę za kogoś.

Wiesz to nie chodzi o odwalenie roboty, ja po prostu tego nie rozumiem, a jak mam płacić, to wolę zapłacić już za warunek, bo tu nawet nie mam pewności, że ktoś mi to napisze, poza tym potrzebuję to na jutro. Aha i to nie jest z matmy, tak w roli ścisłości, z matmą sama dałabym sobie radę. Poza tym tak przypadkiem weszłam na to forum, więc nawet nie wiedziałam, że jak ktoś prosi o pomoc w programie, to już trzeba pisać gdzie indziej, nie mniej jednak dziękuję za "pomoc"

Ech... wpisujesz do Google metoda trapezów, patrzysz na 4 link, patrzysz, a tu... tak! Kod programu! Czy to jest możliwe?
http://edu.i-lo.tarnow.pl/inf/alg/004_int/0003.php

PS. Musisz tylko zamiast ich funkcji x^2+2x podstawić swoją, czyli ln(x), tak? No to po linijce program blablabla... dodajesz linijkę uses Math; i używasz funkcji z tej rozpiski.

Czy da się jeszcze prościej napisać?

Nie wiem wiesz ja tego zupełnie nie rozumiem... Jak dla mnie jest to śmieszne i do niczego niepotrzebne. Poza tym nie myśl, że jestem taka głupia, parę przykładowych programów spisałam, ale każdy jest inny i nie wiem o co w nich chodzi. Mimo wszystko dziekuję za odpowiedź.