Równanie różniczkowe w Pascalu ( plus wykres)

Witam wszystkich.
Przepraszam , że wam zaśmiecam forum ale mam ogromną prośbę. Jako ze nie jestem kumaty w programowaniu zwracam się do was z ogromna prośbą. Otóż nasz kochany profesor ze studiów zadał nam pewien projekt do realizacji. Mamy układ mechaniczny , który mamy zapisać w postaci równania różniczkowego i dla danego równania dokonać obliczeń. Problem w tym ze w projekcie tym ma znaleźć się program komputerowy ( najlepiej pascal) który będzie to równanie liczyć. Poszukałem po necie i trochę poczytałem i napisałem ten program przed chwilką i nie wiem czy jest dobrze. Cichłbym as prosić abyście sprawdzili czy dany program jets prawidłowy ( ewentualnie go zmodyfikowali)...ale najgorsze jest to ze musze także przedstawić obliczenie tego równania na wykresie a nie mam pojęcia jak w Pascalu robi się wykresy. Był bym bardzo wdzięczny gdy jakaś osoba mi pomogła. potrzebuje ten projekt jak najszybciej. Mam już wszystko obliczone i zrobione ale został mi jeszcze właśnie program i wykres.

a oto równanie:

y" = -(Kr/m)y' - Ks/my + Kr/m*x'

gdzie : Kr, Ks, m to dane parametry a x= 1 a w drugim programie x= Asinwt

mam to równanie różniczkowe rozwiązać za pomocą podwójnego kroku całkowania a to mój program:

Program Marcin M;{metoda o podwójnym kroku całkowania}
uses crt;
var
y1,y0,yp0,yp1,zp0,zp1,z0,z1,t,x,Kr,Ks,dt,tk:real;
i,m:integer;
begin

write(?krok cal. Dt=?);readln(dt);
write(?czas obserwacji Tk?);readln(tk);
write(?wydruk wyników m?);readln(m);
writeln;
writeln(?Tk=?,tk4,? ?,?dt=?,dt4);

t:=0;y1:=0;y0:=0;{ustawienie warunków początkowych}
repeat
for i:=0 to m do
begin
x:=1;{obliczenie zmiennych parametrów i funkcji zależnych od czasu}
z0:=y1+Kr/m;z1:=-(Kr/m)y1-(Ks/m)y0;{oblicz wartości na wejściach bloków całkujących}
yp0:=yo+z0dt;yp1:=y1+z1*dt;{oblicz wartości na wyjściach bloków całkujących}

zp0:=yp1+Kr/m,zp1:=-(Kr/m)yp1-(Ks/m)yp0;
y0:=y0+0.5dt(z0+zp0);y1:=y1+0.5dt(z1+zp1);
t:=t+dt;{aktualizacja czasu}
end;
writeln(t6,? ?,y06);{wydruk wyników obliczeń}
until t>=tk;
end.

i drugi dla x=Asinwt

Program Marcin Mann;{metoda o podwójnym kroku całkowania}
uses crt;
var
y1,y0,yp0,yp1,zp0,zp1,z0,z1,w,t,x,Kr,Ks,dt,tk:real;
i,m:integer;
begin

write(?krok cal. Dt=?);readln(dt);
write(?czas obserwacji Tk?);readln(tk);
write(?wydruk wyników m?);readln(m);
writeln;
writeln(?Tk=?,tk4,? ?,?dt=?,dt4);

t:=0;y1:=0;y0:=0;{ustawienie warunków początkowych}
repeat
for i:=0 to m do
begin
x:=(A/w)cos(wt);{obliczenie zmiennych parametrów i funkcji zależnych od czasu}
z0:=y1;z1:=(Kr/m)x-(Kr/m)y1-(Ks/m)y0;{oblicz wartości na wejściach bloków całkujących}
yp0:=yo+z0dt;yp1:=y1+z1dt;{oblicz wartości na wyjściach bloków całkujących}

zp0:=yp1,zp1:=(Kr/m)x(t+dt)-(Kr/m)yp1-(Ks/m)yp0;
y0:=y0+0.5dt(z0+zp0);y1:=y1+0.5dt(z1+zp1);
t:=t+dt;{aktualizacja czasu}
end;
writeln(t6,? ?,y06);{wydruk wyników obliczeń}
until t>=tk;
end.

Wykres w paskalu w trybie dos, więc musisz unit graph użyć.
są tam funkcje initgraphics, lineto, linerel, setcolor(aż 16 sztuk), itp. -> F1 ci pomoże.

Ten pierwszy program nic nie liczy. Nie podałeś wartości Ks,Kr.
Jeśli możesz to opisz krótko co to za zagadnienie mechaniczne.
Na pierwszy rzut oka to jakieś drgania z tłumieniem i wymuszeniem?
PS
Wyniki wyświetlane to narastający czas i same zera.
Jeśli w równaniu występuje x' jako pochodna dx/dt to jeśli x=1 to d(1)/dt=0;
Dużo znaków zapytania!
Czy to jest taka postć ?
d2y/dt2=(-Kr/m)*dy/dt-(Ks/m)*y-(Kr/m)*dx/dt

Witam to jest układ mechaniczny w sumie profesor wymyślił go z głowy to jest jakaś masa, połączona z tłumikiem i sprężyna a równanie które napisałem opisuje ten ruch. Ale równanie jest dobre .Kr i Ks to stałe które mam przyjąć dowolnie ja wiozłem Kr = 40 , Ks = 30 m = 10 wtedy delta tego równania ma przyzwoite pierwiastki r1 = -1 r2 = -3 .Nie wiem co program liczy nie sprawdzałem go jeszcze . W równaniu moje x? mam przyjąć w pierwszym przypadku x=1 .x? jest impulsem , tym się nie przejmuj zrobiłem na kartce obliczenia i schemat blokowy układu całkującego. I równania wyszły jak wyszły prawdopodobnie są dobre. Ja tylko chciałbym byście sprawdzili czy program działa poprawnie czy on liczy w ogóle to a ma liczyć te równania co są podane w programie czyli obliczać wejścia z0 z1 a następnie wyjścia y0 i y1 a następnie ma wykonać krok całkujący czyli zp0 ,zp1 itp...

P.s StPan jeśli masz GG to podaj napisze Ci konkretnie na GG o co biega i podeślę Ci skany tego zadanka i schemat blokowy tego programu oraz algorytm

Zapomniałem podać mój nr GG 5495587

Jeśli temat jest jeszcze aktualny to prześlij na mój adres
[email protected] swój e-mail.
GG nie używam.

Witam...tak temat nadal aktualny z tym ze program już sam poprawiłem i działa teraz kwestia tego jak zrobić by narysował mi pascal wykres dla moich 10 wyników które program mi wylicza. ???StPan wolałbym pisać na GG bo częstotliwość maila nie jest zbyt szybka :p a zależy mi na czasie bo musze mieć to do jutra wieczorem gotowe ....pozdrawiam

To jest program bardziej rozbudowany.
obrazuje drgania wymuszone z tłumieniem tak jak w twoim modelu.
pokazałem jak można bawic się przy różnych danych wejścia.
Jest też częśc graficzna i możliwosc drukowania wynikow.
Do rozwiązania r.różniczkowego użyłem metody Rungego-Kutty.

program oscylator;
uses Crt,Graph;
const n=200;
var 
     mi,tryb                  :Integer;
     p                        :Char;
     m,A,omega,Kr,Ks,
     x0,t0,tk,p0,px,t,x       :Real;

     tb,xb,px1                :array[0..n+1] of Real;
     wydruk    :Text;

function g(t,x,px:Real)             :Real;
   begin
    g:=-(Kr/m)*px-(Ks/m)*x+(A/m)*sin(omega*t);;
   end;
function f(t,x,px:Real)             :Real;
   begin
    f:=1;
   end;

procedure RK2rz2;
  var    h,k1,k2,k3,k4,k,l                :Real;
         s                                :Integer;
begin
xb[0]:=x0;
tb[0]:=t0;
px1[0]:=p0;
h:=tk/n;
  for s:=0 to n do
   begin
    mi:=s+1;
    tb[s]:=tb[0]+s*h;
    t:=tb[s];
    x:=xb[s];
    px:=px1[s];
    k1:=0.5*Sqr(h)*g(t,x,px);
    t:=tb[s]+h/2;
    x:=xb[s]+px1[s]*h/2+k1/4;
    px:=px1[s]+k1/h;
    k2:=0.5*Sqr(h)*g(t,x,px);
    px:=px1[s]+k2/h;
    k3:=0.5*Sqr(h)*g(t,x,px);
    t:=tb[s]+h;
    x:=xb[s]+px1[s]*h+k3;
    px:=px1[s]+2*k3/h;
    k4:=0.5*Sqr(h)*g(t,x,px);
    k:=(1/3)*(k1+k2+k3);
    l:=(1/3)*(k1+2*k2+2*k3+k4);
    tb[mi]:=tb[0]+mi*h;
    px1[mi]:=px1[s]+l/h;
    xb[mi]:=xb[s]+h*px1[s]+k;
   end;
 end;
procedure RK2rz1;
  var    h,k1,k2,k3,k4,k,l                :Real;
         s                                :Integer;
begin
xb[0]:=x0;
tb[0]:=t0;
px1[0]:=p0;
h:=tk/n;
  for s:=0 to n do
   begin
    mi:=s+1;
    tb[s]:=tb[0]+s*h;
    t:=tb[s];
    x:=xb[s];
    px:=px1[s];
    k1:=0.5*Sqr(h)*f(t,x,px);
    t:=tb[s]+h/2;
    x:=xb[s]+px1[s]*h/2+k1/4;
    px:=px1[s]+k1/h;
    k2:=0.5*Sqr(h)*f(t,x,px);
    px:=px1[s]+k2/h;
    k3:=0.5*Sqr(h)*f(t,x,px);
    t:=tb[s]+h;
    x:=xb[s]+px1[s]*h+k3;
    px:=px1[s]+2*k3/h;
    k4:=0.5*Sqr(h)*f(t,x,px);
    k:=(1/3)*(k1+k2+k3);
    l:=(1/3)*(k1+2*k2+2*k3+k4);
    tb[mi]:=tb[0]+mi*h;
    px1[mi]:=px1[s]+l/h;
    xb[mi]:=xb[s]+h*px1[s]+k;
   end;
 end;
procedure model;
begin
Line(20,10,40,10);
Line(30,10,30,20);
MoveTo(30,20);
LineTo(40,30);
LineTo(20,40);
LineTo(40,50);
LineTo(20,60);
LineTo(40,70);
LineTo(30,80);
LineTo(30,90);
Circle(30,100,10);
Line(30,100,30,130);
MoveTo(30,130);
LineTo(35,125);
LineTo(25,125);
LineTo(30,130);
Line(40,100,50,100);
MoveTo(50,100);
LineTo(50,120);
Line(47,120,53,120);
Line(45,110,45,130);
MoveTo(45,130);
LineTo(55,130);
LineTo(55,110);
Line(50,100,80,100);
MoveTo(80,100);
LineTo(80,130); LineTo(85,125); LineTo(75,125);
LineTo(80,130);
Line(50,80,50,60);
MoveTo(50,60);
LineTo(47,65); LineTo(53,65); LineTo(50,60);
Line(30,130,30,180);
OutTextXY(10,135,'mg');
OutTextXY(10,180,'y'); 
Line(50,160,50,135);
MoveTo(50,135);
LineTo(47,140); LineTo(53,140); LineTo(50,135);
OutTextXY(60,150,'Kr*dy/dt');
OutTextXY(60,70,'Ks*y');
OutTextXY(90,120,'F(t)-sila wymuszajaca');
OutTextXY(180,10,'Model Fizyczny Drgan');
OutTextXY(180,230,'Model Matematyczny');
OutTextXY(10,270,'d(dy/dt)/dt=-(Kr/m)*dy/dt-(Ks/m)*y+F(t)/m');
OutTextXY(10,460,'Wcisnij ENTER');
end;
procedure osie;
 var   i,di      :Integer;
begin
Line(0,240,640,240);
Line(0,0,0,480);

 i:=0;
 di:=40;
  while i<=640 do
   begin
    i:=i+di;
    Line(i,240,i,243);
   end;
 i:=0;
 di:=20;
  while i<=480 do
   begin
    i:=i+di;
    Line(0,i,3,i);
   end;
OutTextXY(40,245,'1.0');
OutTextXY(630,245,'t');
OutTextXY(5,220,'0.005');
OutTextXY(5,10,'y');
OutTextXY(10,460,'Wcisnij ENTER');
end;
procedure inekr;
var
    sterownik,blad  :Integer;
begin
sterownik:=Detect;
InitGraph(sterownik,tryb,'C:\Pascal\TP\BGI');{podaj wlasna sciezke do BGI}
blad:=GraphResult;
 if blad<>0
  then
   begin
    Writeln('blad podczas inicjowania trybu graficznego');
    Writeln('kod bladu:',blad);
    Halt;
   end;
SetViewPort(0,0,639,479,ClipOff);
ClearViewPort;
SetBkColor(White);
SetColor(Blue);
SetTextStyle(SmallFont,HorizDir,5);
SetLineStyle(SolidLn,0,NormWidth);
end;
begin
inekr;
model;
Readln;
TextMode(CO80);
TextBackground(Black);
TextColor(White);
Writeln('Opcje impulsu');
Writeln('1. x=1');
Writeln('2. x=A*sin(omega*t)');
 while KeyPressed do p:=ReadKey;
  p:=ReadKey;
 if p='2' then
  begin
   Writeln('podaj mase 1<=m<=10');
   Readln(m);
   Writeln('podaj wsp. sprezystosci 4<=Ks<=40');
   Readln(Ks);
   Writeln('podaj wsp. tlumienia 4<=Kr<=30');
   Readln(Kr);
   Writeln('podaj czas obserwacji tk<=15 s');
   Readln(tk);
   Writeln('podaj amplitude sily wymuszaj. 0.5<=A<=1');
   Readln(A);
   Writeln('podaj czestosc sily wymusz. 1<=omega<=5');
   Readln(omega);
   SetGraphMode(tryb);
   ClearViewPort;
   SetBkColor(White);
   SetColor(Blue);
   SetTextStyle(SmallFont,HorizDir,5);
   SetLineStyle(SolidLn,0,NormWidth);
   osie;

   x0:=0;
   p0:=0;
   t0:=0;
   RK2rz2;
   SetColor(Red);
   MoveTo(0,240);
    for mi:=1 to n do
     begin
      LineTo(Round(40*tb[mi]),Round(240-4000*xb[mi]))
     end;
   Readln;
  end
 else if p='1'then
  begin
   Writeln('podaj mase 1<=m<=10');
   Readln(m);
   Writeln('podaj wsp. sprezystosci 4<=Ks<=40');
   Readln(Ks);
   Writeln('podaj wsp. tlumienia 4<=Kr<=30');
   Readln(Kr);
   Writeln('podaj czas obserwacji tk<=15 s');
   Readln(tk);
   SetGraphMode(tryb);
   ClearViewPort;
   SetBkColor(White);
   SetColor(Blue);
   SetTextStyle(SmallFont,HorizDir,5);
   SetLineStyle(SolidLn,0,NormWidth);
   osie;

   x0:=0;
   p0:=0;
   t0:=0;
   RK2rz1;
   SetColor(Red);
   OutTextXY(50,220,'*1000');
   MoveTo(0,240);
    for mi:=1 to n do
     begin
      LineTo(Round(40*tb[mi]),Round(240-4*xb[mi]))
     end;
   Readln;
  end;
RestoreCrtMode;

 for mi:=1 to n do
  begin
   Delay(30000);
   Write('t= ',tb[mi]:2:3);Write(' ');
   Writeln('y=',xb[mi]:2:4);
  end;
 Writeln('czy drukowac wyniki T/N?');
   while KeyPressed do p:=ReadKey;
   p:=ReadKey;
    if p='t' then
      begin
        Assign(wydruk,'LPT1');
        Rewrite(wydruk);
          for mi:=0 to n do
           begin
             Write(wydruk,'  ','t=',tb[mi]:2:4,'  ');
             Writeln(wydruk,'y=',xb[mi]:2:4);
           end;
       end
     else if p='n' then 

CloseGraph;
end.

Dzięki za fatygę StPan na pewno program sie przyda ale już nie na egzamin właśnie przed chwilą wróciłem zaliczyłem na 4,5 :d ale wielkie dzięki program robi wrażenie :)