Problem mit C programmierung < Numerik < Hochschule < Mathe < Vorhilfe
|
| Status: |
(Frage) beantwortet  | | Datum: | 18:01 So 09.10.2005 | | Autor: | alphaone |
Ich habe ein Problem bei einer C Programmierung. Die Aufgabe ist es die PDG
[mm] u_{t}+u*u_{x}=0 [/mm]
numerisch mit Hilfe einer spektralen Methode für verschiedene Zeitpunkte zu lösen. Mein Programmcode lautet:
| 1: |
| | 2: | #include <catam.h>
| | 3: |
| | 4: | #define N 128
| | 5: | #define Pi 3.1415926535897932384626433832
| | 6: |
| | 7: | double scale(double U_N[N][2]){
| | 8: |
| | 9: | int i;
| | 10: |
| | 11: | for(i=0;i<N;i++){
| | 12: | U_N[i][0]=U_N[i][0]/(double)N;
| | 13: | U_N[i][1]=U_N[i][1]/(double)N;
| | 14: | }
| | 15: | }
| | 16: |
| | 17: | double reposition(double U_N[N][2]){
| | 18: | int i;
| | 19: | double u[N][2];
| | 20: | for(i=0;i<N/2;i++){
| | 21: | u[i][0]=U_N[i][0];
| | 22: | u[i][1]=U_N[i][1];
| | 23: | U_N[i][0]=U_N[i+N/2][0];
| | 24: | U_N[i][1]=U_N[i+N/2][1];
| | 25: | U_N[i+N/2][0]=u[i][0];
| | 26: | U_N[i+N/2][1]=u[i][1];
| | 27: | }
| | 28: | }
| | 29: |
| | 30: | double spectral(double U_N[N][2],double W[N][2]){
| | 31: | int k,j,i;
| | 32: | double X[N][2],u[N][2],u_x[N][2];
| | 33: | for(k=0;k<=N-1;k++){
| | 34: | u[k][0]=U_N[k][0];
| | 35: | u[k][1]=U_N[k][1];
| | 36: | X[k][0]=-(k-N/2)*U_N[k][1];
| | 37: | X[k][1]=(k-N/2)*U_N[k][0];
| | 38: | }
| | 39: | FftCL(u,7,1.0);
| | 40: | FftCL(X,7,1.0);
| | 41: | for(j=0;j<N;j++){
| | 42: | W[j][0]=u[j][0]*X[j][0]-u[j][1]*X[j][1];
| | 43: | W[j][1]=u[j][1]*X[j][0]+u[j][0]*X[j][1];
| | 44: | }
| | 45: | FftCL(W,7,-1.0);
| | 46: | scale(W);
| | 47: | reposition(W);
| | 48: | }
| | 49: |
| | 50: | double f(double x){
| | 51: | double i;
| | 52: | i=-x;
| | 53: | return i;
| | 54: | }
| | 55: |
| | 56: | double Runge(double W[N][2],int p,double dt,double U_N[N][2]){
| | 57: | double a,kU_1[N][2],kU_2[N][2],kU_3[N][2],kU_4[N][2],u[N][2];
| | 58: | int i,j,y;
| | 59: | for(i=1;i<=p+1;i=i+1){
| | 60: | spectral(U_N,W);
| | 61: | for(j=0;j<N;j++){
| | 62: | kU_1[j][0]=dt*f(W[j][0]);
| | 63: | kU_1[j][1]=dt*f(W[j][1]);
| | 64: | u[j][0]=U_N[j][0]+0.5*kU_1[j][0];
| | 65: | u[j][1]=U_N[j][1]+0.5*kU_1[j][1];
| | 66: | }
| | 67: | spectral(u,W);
| | 68: | for(j=0;j<N;j++){
| | 69: | kU_2[j][0]=dt*f(W[j][0]);
| | 70: | kU_2[j][1]=dt*f(W[j][1]);
| | 71: | u[j][0]=U_N[j][0]+0.5*kU_2[j][0];
| | 72: | u[j][1]=U_N[j][1]+0.5*kU_2[j][1];
| | 73: | }
| | 74: | spectral(u,W);
| | 75: | for(j=0;j<N;j++){
| | 76: | kU_3[j][0]=dt*f(W[j][0]);
| | 77: | kU_3[j][1]=dt*f(W[j][1]);
| | 78: | u[j][0]=U_N[j][0]+kU_3[j][0];
| | 79: | u[j][1]=U_N[j][1]+kU_3[j][1];
| | 80: | }
| | 81: | spectral(u,W);
| | 82: | for(j=0;j<N;j++){
| | 83: | kU_4[j][0]=dt*f(W[j][0]);
| | 84: | kU_4[j][1]=dt*f(W[j][1]);
| | 85: | U_N[j][0]=(U_N[j][0]+(double)1/6*(kU_1[j][0]+2*kU_2[j][0]+2*kU_3[j][0]+kU_4[j][0]));
| | 86: | U_N[j][1]=(U_N[j][1]+(double)1/6*(kU_1[j][1]+2*kU_2[j][1]+2*kU_3[j][1]+kU_4[j][1]));
| | 87: | }
| | 88: | }
| | 89: | }
| | 90: |
| | 91: | double load(double U_N[N][2]){
| | 92: | int i,l;
| | 93: | l=1;
| | 94: | for(i=0;i<N/4;i++){
| | 95: | U_N[i][0]=sin(2*Pi*i/(double)N);
| | 96: | U_N[i][1]=0;
| | 97: | }
| | 98: | U_N[N/4][0]=1;
| | 99: | U_N[N/4][1]=0;
| | 100: | for(i=N/4+1;i<=N/2;i++){
| | 101: | U_N[i][0]=U_N[N/4-l][0];
| | 102: | U_N[i][1]=0;
| | 103: | l=l+1;
| | 104: | }
| | 105: | for(i=N/2+1;i<=3*N/4;i++){
| | 106: | U_N[i][0]=-U_N[i-N/2][0];
| | 107: | U_N[i][1]=0;
| | 108: | }
| | 109: | for(i=3*N/4+1;i<N;i++){
| | 110: | U_N[i][0]=-U_N[i-N/2][0];
| | 111: | U_N[i][1]=0;
| | 112: | }
| | 113: | }
| | 114: |
| | 115: |
| | 116: | int MainCL(void){
| | 117: |
| | 118: | double U_N[N][2],a,q,W[N][2],t,dt,u[N],x[N],y[100000],z[100000];
| | 119: | int i,p,l,e;
| | 120: | p=100;
| | 121: | printf("Please enter the time point\n");
| | 122: | scanf("%lf",&t);
| | 123: |
| | 124: | dt=t/(double)p;
| | 125: | load(U_N);//initial distrbution: U(x)=sin(x) at t=0//
| | 126: | FftCL(U_N,7,-1.0);//fourier transform//
| | 127: | scale(U_N);//scaling of fourier transform not done above- so done now//
| | 128: | reposition(U_N);//exchange lower half of wave numbers with upper half//
| | 129: | Runge(W,p,dt,U_N);//perform runge kutta integration//
| | 130: | FftCL(U_N,7,1.0);//inverse fourier//
| | 131: |
| | 132: | for(i=0;i<N;i++){
| | 133: | x[i]=2*Pi*i/(double)N;
| | 134: | u[i]=U_N[i][0];
| | 135: | }
| | 136: | for(i=0;i<1000;i++){
| | 137: | y[i]=sin(i/(float)100);
| | 138: | z[i]=i/(float)100+t*sin(i/(float)100);
| | 139: | }
| | 140: |
| | 141: | XYCurveCL(x,u,N,1,DOT,RedCC,AUTOAXES);
| | 142: | XYCurveCL(z,y,1000,1,DOT,BlueCC,PRESET);
| | 143: | }
|
wobei XYCurveCL eine vorgegebene Funktion die x gegen u plottet und FftCL eine Funktion ist die einen Fast Fourier Transform durchführt, beide Funktionen funktionieren soweit mir bekannt ist korrekt. y[i] und z[i] stellen die exakte Lösung parametrisiert dar und [mm] U_N[i][0] [/mm] gibt die Lösung, die mit hilfe von spektraler Methode und Runge Kutta berechnet wurde dar. Das Problem ist nun, dass die Form der Lösung zwar sehr nah an der exakten Lösung ist, jeder zweite wert in der numerischen Lösung aber das falsche Vorzeichen hat. Selbstverständlich könnte ich jetzt jeden zweiten Term mit -1 multiplizieren, jedoch wäre es sicher zu bevorzugen zunächst einmal zu verstehen wo der Vorzeichendreher herkommt, was mir bis jetzt aber nicht geglückt ist. Sollte irgendjemand eine Idee haben, wäre ich sehr dankbar.
|
|
| |
|
Hallo alphaone,
ich kann deinem Programm zwar nicht entnehmen, wie es im einzelnen arbeitet, aber mich hat die Funktion
f(double x),
die ja nur das Vorzeichen vertauscht, etwas irritiert. Vielleicht kannst du sie aus dem Programm entfernen und findest dabei den Fehler. Prinzipiell würde ich versuchen, alle Casts zu vereinheitlichen, d.h. entweder alles auf float oder alles auf double. Das ist aber wahrscheinlich nicht das eigentliche Problem.
Meine Fragen sind vielmehr:
1. ist der Vorzeichenwechsel evtl. gar kein Fehler?
2. kannst du das Programm etwas ausführlicher kommentieren?
3. kannst du den Code so formatieren, dass er auch im Forum gut lesbar
ist?
Was passiert, wenn du die Schrittweite deiner Diskretisierung halbierst/verdoppelst? Tritt dann immer noch der VZW bei jeder zweiten Zahl auf?
Hugo
|
|
|
|
