Jak zaimplementować pochodną funkcji
Oto mój pomysł
1. Utworzyć trzy ciągi korzystając z ilorazów różnicowych
a) \(\displaystyle{ \frac{f \left( x+h\right)-f \left( x\right) }{h}}\)
Iloraz różnicowy przedni
b) \(\displaystyle{ \frac{f \left( x+h\right)-f \left( x-h\right) }{2h}}\)
Iloraz różnicowy centralny
c) \(\displaystyle{ \frac{f \left( x\right)-f \left( x-h\right) }{h}}\)
Iloraz różnicowy wsteczny
2. Obliczyć granice tych ciągów i jeżeli te granice są równe
to granica ta jest równa pochodnej funkcji w punkcie
I jak to zaimplementować w językach takich jak Pascal czy C
i czy mój pomysł jest dobry
Różniczkowanie numeryczne
-
spajder
- Użytkownik
- Posty: 735
- Rejestracja: 7 lis 2005, o 23:56
- Płeć: Mężczyzna
- Lokalizacja: Łódź
- Podziękował: 2 razy
- Pomógł: 133 razy
Różniczkowanie numeryczne
W ten sposób nic nie zdziałasz, bo nie jesteś w stanie numerycznie wyliczyć granicy. Poszukaj:
1. metoda Eulera (prosta)
2. metoda Rungego-Kutty (mnie ciężko było załapać)
jest oczywiście jeszcze sporo innych, ale te dwie są popularne i w miarę łatwo się je implementuje
1. metoda Eulera (prosta)
2. metoda Rungego-Kutty (mnie ciężko było załapać)
jest oczywiście jeszcze sporo innych, ale te dwie są popularne i w miarę łatwo się je implementuje
-
Mariusz M
- Użytkownik
- Posty: 6909
- Rejestracja: 25 wrz 2007, o 01:03
- Płeć: Mężczyzna
- Lokalizacja: 53°02'N 18°35'E
- Podziękował: 2 razy
- Pomógł: 1246 razy
Różniczkowanie numeryczne
Spajder te metody które podałeś dotyczą numerycznego rozwiązywania równań
różniczkowych
Mógłbyś mi objaśnić poniższy kod podobno tutaj biorą granicę z ilorazu centralnego
Ty napisałeś że to niemożliwe
różniczkowych
Mógłbyś mi objaśnić poniższy kod podobno tutaj biorą granicę z ilorazu centralnego
Ty napisałeś że to niemożliwe
Kod: Zaznacz cały
program NUMERICALDIFFERENTIATIONLIMIT;
{--------------------------------------------------------------------}
{ Alg6'12.pas Pascal program for implementing Algorithm 6.1-2 }
{ }
{ NUMERICAL METHODS: Pascal Programs, (c) John H. Mathews 1995 }
{ To accompany the text: }
{ NUMERICAL METHODS for Math., Science & Engineering, 2nd Ed, 1992 }
{ Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. }
{ Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6 }
{ Prentice Hall, International Editions: ISBN 0-13-625047-5 }
{ This free software is compliments of the author. }
{ E-mail address: in%"mathews@fullerton.edu" }
{ }
{ Algorithm 6.2 (Differentiation Using Extrapolation). }
{ Section 6.1, Approximating the Derivative, Page 327 }
{--------------------------------------------------------------------}
uses
crt;
const
Delta = 1E-7;
Tol = 1E-7;
Small = 1E-7;
Big = 1E6;
Max = 32;
MaxN = 15;
MaxV = 34;
FunMax = 9;
type
VECTOR = array[0..MaxV] of real;
MATRIX = array[0..MaxN, 0..MaxN] of real;
LETTER = string[1];
Status = (Done, Working);
DoSome = (Go, More, Stop);
LETTERS = string[200];
var
FunType, Inum, J, K, M, Meth, N, NS, Sub: integer;
DL, DR, DS, Error, ErrorS, H, H0, HL, HR, RelErr, Rnum, X: real;
D, DD, E, F1, F2, HH: VECTOR;
Ans: CHAR;
Stat: Status;
DoMo: DoSome;
Mess: LETTERS;
function F (X: real): real;
begin
case FunType of
1:
F := SIN(X);
2:
F := COS(X);
3:
F := EXP(X);
4:
begin
if X = 0 then
F := 1E-38
else
F := LN(ABS(X));
end;
5:
F := ABS(X);
6:
F := SQRT(ABS(X));
7:
begin
if X = 0 then
F := 0
else
begin
if X > 0 then
F := EXP(LN(X) / 3)
else
F := -EXP(LN(ABS(X)) / 3);
end;
end;
8:
begin
if X = 0 then
F := 0
else
begin
if X > 0 then
F := EXP(LN(X * X) / 3)
else
F := EXP(LN(ABS(X * X)) / 3);
end;
end;
9:
F := ABS(X) - X;
end;
end;
procedure PRINTFUNCTION (FunType: integer);
begin
case FunType of
1:
WRITELN('F(X) = SIN(X)');
2:
WRITELN('F(X) = COS(X)');
3:
WRITELN('F(X) = EXP(X)');
4:
WRITELN('F(X) = LN(X)');
5:
WRITELN('F(X) = |X|');
6:
WRITELN('F(X) = |X|^1/2');
7:
WRITELN('F(X) = X^1/3');
8:
WRITELN('F(X) = X^2/3');
9:
WRITELN('F(X) = |X| - X');
end;
end;
procedure TAKELIMIT ( {FUNCTION F(X:real): real;}
var D, E: VECTOR; var X, H: real; var M: integer);
var
J, Min, N: integer;
R: VECTOR;
begin
H := H0;
Min := 2;
D[0] := 0.5 * (F(X + H) - F(X - H)) / H;
F2[0] := F(X + H); {Storing for details}
F1[0] := F(X - H);
HH[0] := H;
for J := 1 to 2 do
begin
H := H / 2;
HH[J] := H;
D[J] := 0.5 * (F(X + H) - F(X - H)) / H;
F2[J] := F(X + H); {Storing for details}
F1[J] := F(X - H);
E[J] := ABS(D[J] - D[J - 1]);
R[J] := 2 * E[J] / (ABS(D[J]) + ABS(D[J - 1]) + Small);
if (ABS(D[J]) < 0.001) then
Min := 8;
end;
N := 1;
while (((E[N] > E[N + 1]) or (R[N] > Tol)) and (N < Max) and (ABS(D[N]) < Big)) or (N < Min) do
begin
H := H / 2;
HH[N + 2] := H;
D[N + 2] := 0.5 * (F(X + H) - F(X - H)) / H;
F2[N + 2] := F(X + H); {Storing for details}
F1[N + 2] := F(X - H);
E[N + 2] := ABS(D[N + 2] - D[N + 1]);
R[N + 2] := 2 * E[N + 2] / (ABS(D[N + 2]) + ABS(D[N + 1]) + Small);
if (ABS(D[N + 2]) < 0.001) then
Min := 8;
N := N + 1;
end;
H := 2 * H;
M := N;
end;
procedure EXTRAPOLATE (var X, DE, HL, HR, Error: real; var DD: VECTOR; var N: integer);
const
MaxN = 15;
type
MATRIX = array[0..MaxN, 0..MaxN] of real;
var
J, K: integer;
RelErr: real;
D: MATRIX;
function E4 (L: integer): real;
var
I: integer;
P: real;
begin
P := 1;
for I := 1 to L do
P := P * 4;
E4 := P;
end;
begin
J := 1;
Error := 1;
RelErr := 1;
D[0, 0] := (F(X + HR) - F(X - HL)) / (HL + HR);
DD[0] := D[0, 0];
HH[0] := HL + HR;
while (RelErr > Tol) and (Error > Delta) and (J <= MaxN) do
begin
HL := HL / 2;
HR := HR / 2;
HH[J] := HL + HR;
D[J, 0] := (F(X + HR) - F(X - HL)) / (HL + HR);
for K := 1 to J do
D[J, K] := D[J, K - 1] + (D[J, K - 1] - D[J - 1, K - 1]) / (E4(K) - 1);
Error := ABS(D[J, J] - D[J - 1, J - 1]);
RelErr := 2 * Error / (ABS(D[J, J]) + ABS(D[J - 1, J - 1]) + Small);
DD[J] := D[J, J];
N := J;
J := J + 1;
end;
DE := D[N, N];
end;
procedure GETFUN (var FunType: integer);
var
K: integer;
begin
CLRSCR;
WRITELN(' Choose your function:');
WRITELN;
for K := 1 to FunMax do
begin
WRITE(' <', K : 2, ' > ');
PRINTFUNCTION(K);
WRITELN;
end;
WRITELN;
WRITE(' Select < 1 - ', FunMax : 1, ' > ? ');
FunType := 1;
READLN(FunType);
if FunType < 1 then
FunType := 1;
if FunType > FunMax then
FunType := FunMax;
end;
procedure GETX (Meth: integer; var X: real);
begin
CLRSCR;
WRITELN;
case Meth of
1:
begin
WRITELN(' The derivative will be approximated using the limit of a');
WRITELN;
WRITELN('sequence of symmetric difference quotients:');
WRITELN;
WRITELN(' -n ');
WRITELN(' lim D = lim [f(x+h ) - f(x-h )]/2h where h = 2 .');
WRITELN('n->oo n n->oo n n n n ');
WRITELN;
WRITELN;
WRITELN('Computation continues until |D - D | => |D - D | or until');
WRITELN(' N+1 N N N-1');
WRITELN;
WRITELN('consecutive approximations agree to at least six decimal places.');
end;
2:
begin
WRITELN('Richardson`s extrapolation is applied to the sequence of difference quotients');
WRITELN;
WRITELN(' -k ');
WRITELN(' { [f(x+h ) - f(x-h )]/2h : where h = 2 }');
WRITELN(' k k k k ');
WRITELN;
WRITELN('and an approximation to the derivative of the function is determined.');
WRITELN;
WRITELN;
WRITELN('Extrapolation continues until consecutive approximations');
WRITELN;
WRITELN;
WRITELN('agree to at least six decimal places.');
end;
end;
WRITELN;
WRITELN(' You chose to find the derivative of the function');
WRITELN;
WRITE(' ');
PRINTFUNCTION(FunType);
WRITELN;
Mess := ' ENTER the value X = ';
X := 0;
WRITE(Mess);
READLN(X);
Mess := ' ENTER initial size for H = ';
H0 := 1;
WRITE(Mess);
READLN(H0);
if H0 < 0 then
H0 := 1;
if H0 < 1E-9 then
H0 := 1E-9;
end;
procedure RESULTS (D, E: VECTOR; X, H: real; M: integer);
var
B, DL, DLE, DR, DRE: real;
begin
CLRSCR;
WRITELN;
WRITELN(' The limit of the symmetric difference quotients was found:');
WRITELN;
WRITELN(' f`(x) = lim [ f(x+h) - f(x-h) ] / 2h ,');
WRITELN(' h -> 0 ');
WRITELN;
WRITE(' where ');
PRINTFUNCTION(FunType);
WRITELN;
WRITELN(' After ', M : 2, ' terms were computed, the');
WRITELN;
WRITELN('the step size was reduced to h = ', H : 15 : 7);
WRITELN;
WRITELN(' The approximation for the value of the derivative is:');
WRITELN;
WRITELN(' f`(', X : 15 : 7, ' ) = ', D[M] : 15 : 7);
WRITELN;
WRITELN('The accuracy is +-', E[M] : 15 : 7);
DL := (F(X) - F(X - H)) / H;
DR := (F(X + H) - F(X)) / H;
B := ABS(D[M - 1]) + ABS(DL) + ABS(DR) + 0.01;
DLE := ABS(DL - D[M]) / B;
DRE := ABS(DR - D[M]) / B;
if (DL * DR < 0) then
if (ABS(DL) < 0.005) and (ABS(DR) < 0.005) then
begin
DLE := 0;
DRE := 0;
end;
if (DLE > 0.001) or (DRE > 0.001) then
begin
WRITELN;
WRITELN(' The symmetric derivative does not agree with');
WRITELN(' the left-hand and right-hand derivatives.');
WRITELN;
WRITELN(' Left-hand [ f(x) - f(x-h) ]/h = ', DL : 15 : 7);
WRITELN;
WRITELN(' Right-hand [ f(x+h) - f(x) ]/h = ', DR : 15 : 7);
WRITELN;
WRITE('Perhaps f`(x) does not exist. Press the <ENTER> key.');
READLN(Ans);
WRITELN;
end;
if (ABS(D[M]) > Big) then
begin
WRITELN;
WRITELN(' The difference quotient seems large.');
WRITELN;
WRITELN('Perhaps f`(x) is infinite. Press the <ENTER> key.');
READLN(Ans);
WRITELN;
end;
end;
procedure RESULTE (X, DS, DL, DR, H, HL, HR, Error: real; NS: integer);
var
B: real;
begin
CLRSCR;
WRITELN;
WRITELN('Richardson`s extrapolation was used on the sequence of difference quotients:');
WRITELN;
WRITELN(' -k ');
WRITELN(' { [f(x+h ) - f(x-h )]/2h : where h = 2 }');
WRITELN(' k k k k ');
WRITELN;
WRITELN('and an approximation was found for the derivative of the function:');
WRITELN;
PRINTFUNCTION(FunType);
WRITELN;
WRITELN(' After ', NS : 2, ' iterations, the');
WRITELN;
WRITELN('step size was reduced to h = ', H : 15 : 7);
WRITELN;
WRITELN('An approximation for the value of the derivative is;');
WRITELN;
WRITELN('f`(', X : 15 : 7, ' ) = ', DS : 15 : 7);
WRITELN;
WRITELN('The accuracy is +-', Error : 15 : 7);
B := 1 + ABS(DL) + ABS(DR);
if ABS(DR - DL) > 0.001 * B then
begin
WRITELN;
WRITELN('The left-hand and right-hand extrapolations are not equal.');
WRITELN;
WRITELN(' lim [ f(x) - f(x-h) ]/h = ', DL : 15 : 7);
WRITELN('h -> 0');
WRITELN(' lim [ f(x+h) - f(x) ]/h = ', DR : 15 : 7);
WRITELN('h -> 0');
WRITE('Perhaps f`(x) does not exist. Press the <ENTER> key.');
READLN(Ans);
WRITELN;
end
end;
procedure MESSAGE (var Meth: integer);
begin
CLRSCR;
WRITELN(' TECHNIQUES FOR NUMERICAL DIFFERENTIATION');
WRITELN;
WRITELN;
WRITELN(' Choose a method for finding the numerical derivative.');
WRITELN;
WRITELN;
WRITELN(' < 1 > Use the limit of symmetric difference quotients');
WRITELN;
WRITELN(' to find an approximation to f`(x):');
WRITELN;
WRITELN(' f`(x) = lim [ f(x+h) - f(x-h) ]/2h');
WRITELN(' h -> 0 ');
WRITELN;
WRITELN;
WRITELN(' < 2 > Use Richardson`s extrapolation applied');
WRITELN;
WRITELN(' to the sequence of difference quotients:');
WRITELN;
WRITELN(' -k ');
WRITELN(' { [f(x+h ) - f(x-h )]/2h : where h = 2 }');
WRITELN(' k k k k ');
WRITELN;
WRITELN;
Mess := ' Select the method. < 1 or 2 > ';
Meth := 0;
WRITE(Mess);
READLN(Meth);
if Meth < 1 then
Meth := 1;
if Meth > 2 then
Meth := 2;
end;
procedure PRINTAPPROXS;
var
K: integer;
begin
CLRSCR;
if Meth = 1 then
begin
WRITELN;
WRITELN(' D = [ f(x+h ) - f(x+h ) ] / 2 h ');
WRITELN(' k k k k ');
WRITELN('--------------------------------------------------------------------------------');
for K := 0 to M do
begin
if 0 <= F1[K] then
WRITELN(D[K] : 15 : 7, ' = [', F2[K] : 15 : 7, ' - ', F1[K] : 15 : 7, '] /', 2 * HH[K] : 11 : 7)
else
WRITELN(D[K] : 15 : 7, ' = [', F2[K] : 15 : 7, ' + ', ABS(F1[K]) : 15 : 7, '] /', 2 * HH[K] : 11 : 7);
WRITELN;
if K mod 11 = 9 then
begin
WRITE(' Press the <ENTER> key. ');
READLN(Ans);
WRITELN;
WRITELN;
end;
end;
WRITELN;
WRITE(' Press the <ENTER> key. ');
READLN(Ans);
WRITELN;
WRITELN;
WRITELN;
end;
if Meth = 2 then
begin
WRITELN;
WRITELN(' k h D ');
WRITELN(' k k ');
WRITELN(' ---------------------------------------------------');
WRITELN;
for K := 0 to NS do
begin
WRITELN(' ', K : 2, ' ', HH[K] / 2 : 15 : 7, ' ', DD[K] : 15 : 7);
WRITELN;
if K mod 11 = 9 then
begin
WRITE(' Press the <ENTER> key. ');
READLN(Ans);
WRITELN;
WRITELN;
end;
end;
WRITELN;
WRITE(' Press the <ENTER> key. ');
READLN(Ans);
WRITELN;
WRITELN;
end;
end;
begin {Begin Main Program}
MESSAGE(Meth);
DoMo := Go;
while (DoMo = Go) or (DoMo = More) do
begin
if DoMo = More then
begin
WRITELN;
WRITE('Do you want to change the method ? <Y/N> ');
READLN(Ans);
WRITELN;
if (Ans = 'y') or (Ans = 'Y') then
MESSAGE(Meth);
end;
GETFUN(FunType);
Stat := Working;
while Stat = Working do
begin
GETX(Meth, X);
case Meth of
1:
begin
TAKELIMIT(D, E, X, H, M);
RESULTS(D, E, X, H, M);
end;
2:
begin
H := H0;
HL := H0;
HR := 0;
EXTRAPOLATE(X, DL, HL, HR, Error, DD, N);
HL := 0;
HR := H0;
EXTRAPOLATE(X, DR, HL, HR, Error, DD, N);
HL := H0;
HR := H0;
EXTRAPOLATE(X, DS, HL, HR, ErrorS, DD, NS);
RESULTE(X, DS, DL, DR, H, HL, HR, ErrorS, NS);
end;
end;
WRITELN;
WRITE('Want to see all the approximations ? <Y/N> ');
READLN(Ans);
WRITELN;
if (Ans = 'Y') or (Ans = 'y') then
PRINTAPPROXS;
WRITELN;
WRITE('Want evaluate f`(x) at another point ? <Y/N> ');
READLN(Ans);
WRITELN;
if (Ans <> 'Y') and (Ans <> 'y') then
Stat := Done;
end;
WRITELN;
WRITE('Do you want to try another function ? <Y/N> ');
READLN(Ans);
WRITELN;
if (Ans = 'Y') or (Ans = 'y') then
DoMo := More
else
DoMo := Stop;
end;
end.
{End of Main Program}
-
xiikzodz
- Użytkownik
- Posty: 1874
- Rejestracja: 4 paź 2008, o 02:13
- Płeć: Kobieta
- Lokalizacja: Lost Hope
- Podziękował: 28 razy
- Pomógł: 502 razy
Różniczkowanie numeryczne
Mozna z uzyciem ilorazu roznicowego. Na przyklad pochodna funkcji \(\displaystyle{ f(x)}\) w punkcie \(\displaystyle{ x_0=0}\) mozna przyblizyc wartoscia \(\displaystyle{ \frac{f(x_0+0.01)-f(x)}{0.01}}\) i bedzie to dobre przyblizenie dla \(\displaystyle{ f(x)=x^2}\), lecz marne dla \(\displaystyle{ f(x)=(2x)^{10}}\). Znaczy sie z jednej strony chcemy w mozliwie malej liczbie krokow dobrze przyblizyc pochodna, a z drugiej strony chcemy uniknac bledow numerycznych wynikajacych np. z dzielenia przez liczby bliskie zeru. Ponizej ustalamy sobie zawczasu dokladnosc. Nie jest to taka dokladnosc w sensie bledu w wyznaczeniu pochodnej, tylko nominalna/umowna dokladnosc systemu pozycyjnego w maszynie.
Liczymy f'(x):
Ciekawe, czy dziala...
Liczymy f'(x):
Kod: Zaznacz cały
const dokladnosc=0.01;
float diff(x)
{
delta=1;
aktualny=1;
poprzedni=0;
while (delta>dokladnosc)
{
poprzedni=aktualny;
aktualny=f(x+delta)+aktualny-aktualny*delta;
delta=min(delta/2,abs(aktualny-poprzedni));
}//while
return aktualny;
}//diff