MATLAB adaptives Simpsonverf. < Matlab < Mathe-Software < Mathe < Vorhilfe
|
Status: |
(Frage) beantwortet | Datum: | 22:37 So 06.07.2014 | Autor: | riju |
Hallo,
ich soll das adaptive Simpsonverfahren mit Hilfe von MATLAB programmieren.
Allerdings erhalte ich eine Endlossschleife.
Kann mir jemand helfen und sagen, warum ich die erhalte? Ich finde meinen Fehler einfach nicht.
Hier mein Algorithmus:
function [ T ] = simpson( a,b,tol )
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
x0=a;
x2=(a+b)/2;
x4=b;
f0=fun(x0);
f2=fun(x2);
f4=fun(x4);
n=0;
T=0;
E=0;
h=x4-x0;
hmin=((b-a)*2^(-24));
maxit=5;
it=1;
disp('x0 x4 T E');
while(h>hmin) && (it<maxit)
x1=(x0+x2)/2;
x3=(x2+x4)/2;
f1=fun(x1);
f3=fun(x3);
h=x4-x0
S=h*(f0+4*f2+f4)/6;
dS=h*(f0+4*f1+2*f2+4*f3+f4)/12;
if (abs(dS-S)>tol*h) && (h>hmin)
n=n+1;
K1(n)=x3;
K2(n)=x4;
K3(n)=f3;
K4(n)=f4;
x4=x2;
x2=x1;
f4=f2;
f2=f1;
else
T=T+dS+(dS-S)/15;
E=E+abs(dS-S)/15;
%ausg=[x(1),x(5),T,E];
%disp('%12.6f %12.6f %12.6f %12.6f [mm] \n',ausg);
[/mm]
if n>0
x0=x4;
x2=K1(n);
x4=K2(n);
f0=f4;
f2=K3(n);
f4=K4(n);
n=n-1;
end
end
it=it+1;
end
end
function f=fun(x)
[mm] f=sqrt(x^3+x+1);
[/mm]
%f=log(x)/x;
end
Vielen Dank im Voraus.
Liebe Grüße
riju
|
|
|
|
Hallo riju,
ich habe deinen Code so wie er ist kopiert und kann keine Endlosschleife entdecken. Was für Parameter übergibst du deiner Funktion?
Noch eine kleine Bitte: Code sollte immer in Code-Blöcke [code][/code] eingeschlossen werden.
Grüße
sijuherm
|
|
|
|
|
Status: |
(Frage) beantwortet | Datum: | 11:14 Mo 07.07.2014 | Autor: | riju |
Hi, also so gibt es keine Endlossschleife. Das stimmt. Ich hatte da noch was eingebaut, dass keine Endlosschleife entsteht. Hier nochmal der Code mit Endlosschleife
1: |
| 2: | function [ T ] = simpson( a,b,tol )
| 3: | %UNTITLED2 Summary of this function goes here
| 4: | % Detailed explanation goes here
| 5: |
| 6: | x0=a;
| 7: | x2=(a+b)/2;
| 8: | x4=b;
| 9: | f0=fun(x0);
| 10: | f2=fun(x2);
| 11: | f4=fun(x4);
| 12: | n=0;
| 13: | T=0;
| 14: | E=0;
| 15: | h=x4-x0;
| 16: | hmin=((b-a)*2^(-24));
| 17: |
| 18: |
| 19: |
| 20: | disp('x0 x4 T E');
| 21: |
| 22: | while(h>hmin)
| 23: | x1=(x0+x2)/2;
| 24: | x3=(x2+x4)/2;
| 25: | f1=fun(x1);
| 26: | f3=fun(x3);
| 27: | h=x4-x0
| 28: | S=h*(f0+4*f2+f4)/6;
| 29: | dS=h*(f0+4*f1+2*f2+4*f3+f4)/12;
| 30: | if (abs(dS-S)>tol*h) && (h>hmin)
| 31: | n=n+1;
| 32: | K1(n)=x3;
| 33: | K2(n)=x4;
| 34: | K3(n)=f3;
| 35: | K4(n)=f4;
| 36: | x4=x2;
| 37: | x2=x1;
| 38: | f4=f2;
| 39: | f2=f1;
| 40: | else
| 41: | T=T+dS+(dS-S)/15;
| 42: | E=E+abs(dS-S)/15;
| 43: |
| 44: | %ausg=[x(1),x(5),T,E];
| 45: | %disp('%12.6f %12.6f %12.6f %12.6f $ \n',ausg); $
| 46: | if n>0
| 47: | x0=x4;
| 48: | x2=K1(n);
| 49: | x4=K2(n);
| 50: | f0=f4;
| 51: | f2=K3(n);
| 52: | f4=K4(n);
| 53: | n=n-1;
| 54: | end
| 55: | end
| 56: |
| 57: |
| 58: | end
| 59: | end
| 60: |
| 61: | function f=fun(x)
| 62: | f=sqrt(x^3+x+1);
| 63: | %f=log(x)/x;
| 64: | end
|
Kannst vielleicht bitte jemand nochma ldrüber schauen?
Liebe Grüße
riju
|
|
|
|
|
Also die einfache Erklärung ist, dass dein h nach wenigen Iterationen den Wert nicht mehr verändert und somit immer größer als hmin ist. Aber das hast du sicherlich auch selebr gemerkt, die Frage ist jetzt natürlich, warum das so ist.
h wird in Zeile 27 aktualisiert und damit indirekt in Zeile 36 bzw 47-49. Hier widerum liegt das Hauptproblem, dass die Bedingung in Zeile 30 nach wenigen Iterationen nicht mehr erfüllt wird und damit die Zeile 36 direkt bzw 47-49 indirekt nicht mehr erreicht wird. Die Ursache hierfür widerum ist, dass S und dS zu ähnlich sind (2 Varianten der Simpsonregel?)
Leider kann ich dir nicht zu 100% sagen, ob du nen Fehler im Algorithmus drin hast (dazu müsste ich mich noch intensiver in das Thema Simpsonverfahren einlesen) oder ob "nur" deine Abbruchbedingungen zu eng gesetzt sind.
Als Abhilfe hast du ja schon von selbst die maximale Zahl der Iterationen begrenzt, allerdings solltest du das Ergebnis, was dabei herauskommt, hinterfragen.
|
|
|
|
|
Status: |
(Frage) überfällig | Datum: | 16:48 Mo 07.07.2014 | Autor: | riju |
Die Abbruchbedingung ist mir leider so vorgegeben, genauso wie hmin. Vielleicht hat ja noch jemand anderes eine Idee.
Ich hab den Algorithmus auf Papier. Meiner Ansicht nach, hab ich den auch richtig "abgeschrieben".
Liebe Grüße
|
|
|
|
|
Status: |
(Mitteilung) Reaktion unnötig | Datum: | 17:20 Mi 09.07.2014 | Autor: | matux |
$MATUXTEXT(ueberfaellige_frage)
|
|
|
|