Giải phương trình vi phân chuyển động của cơ hệ

đồ án tin học trong cơ khí” là môn học quan trọng và cần thiết đối với sinh viên ngành cơ tin kỹ thuật, nó giúp cho sinh viên làm quen với công việc lập trình bằng máy tính để tính toán thiết kế trong cơ khí và dễ dàng tiếp cận công việc khi ra trường. Đề tài của em là giải phương trình vi phân chuyển động của cơ hệ gồm các phần: Phần I: Giới thiệu về hệ phương trình vi phân cấp một và hệ phương trình vi phân cấp cao áp dụng cho các bài toán cơ học. Phần II: Mô tả về cấu tạo của con lắc eliptic và thiết lập phương trình vi phân chuyển động của con lắc. Phần III: Là phần viết chương trình để giải quyết bài toán về hệ phương trình vi phân cấp một và cấp hai. Chương trình này gồm có 2 module Molule1: Giải hệ phương trình vi phân cấp 1 bằng phương pháp RUNGE_KUTTA Module2: Giải hệ phương trình vi phân cấp 2 và áp dụng cho cơ hệ con lắc ELIPTIC. Mục đích của đề tài nhằm khảo sát các chuyển động của 1 cơ hệ nào đó về vị trí, vận tốc…Qua đó vẽ đồ thị mô tả sự dao động của cơ hệ và đưa ra nhận xét. Phần thuyết minh I/ Giải gần đúng hệ phương trình vi phân thường bằng phương pháp RUNGE_KUTTA Xét hệ gồm n phương trình vi phân cấp 1 y’1= f1(x,y1,y2,……,yn) y’2= f2(x,y1,y2,……,yn) (1) ……………………… y’n= fn(x,y1,y2,……,yn) Với y1, y2,……, yn đều là các hàm số của x: yi=yi(x) i =1,2,….,n Ta có thể viết gọn lại như sau: y’i= fi(x,y1,y2,……,yn) Với các điều kiện đầu: yi(x0)=yi0 i =1,2,……,n Ta có thể viết dưới dạng ma trận: y’= f(x,y) Với các điều kiện đầu: y(x0) = y0 Trong đó: y = [y1,y2,……,yn]T y’= [y’1, y’2,…., y’n]T f = [f1, f2,……., fn]T Để giải hệ trên ta có nhiều phương pháp tuy nhiên ở đây ta sử dụng phương pháp RUNGE_KUTTA Biểu thức hồi quy cho phương pháp RUNGE_KUTTA r điểm được viết như sau: zi+1 = z1 + h(pr1k1 + pr2k2 +…….+ prrkr ) Với k1 = f(xi,zi) k2 = f(xi + a2h, zi + hb21k1) ……… kr = f[xi + ar h, zi + h(br1k1 + br2k2 + …+ br,r-1kr-1)] Trong đó prm, m = 1,2…..,r và am, bmj, m = 2,3,…..,r và j = 1,2,……., r-1 là các hằng số được xác định tuỳ theo từng phương pháp. Phương pháp RUNGE_KUTTA là phương pháp một nút vì để tính giá trị của nút sau chỉ cần tính giá trị của nút trước đó. Phương pháp hay dùng nhất là phương pháp một nút 4 điểm với độ chính xác cao vì quá trình tính toán tương đối đơn giản. Sơ đồ của phương pháp RUNGE_KUTTA một nút 4 điểm được viết như sau: yi+1= yi + h(k1+ 2k2 + 2k3 + k4)/6 Với k1= f(xi,yi) k2= f(xi+h/2,yi+hk1/2) k3= f(xi+h/2,yi+hk2/2) k4= f(xi+h,yi+hk3) Trong các vấn đề kỹ thuật ta còn gặp bài toán ở dạng phương trình vi phân cấp cao: y(n)= f(x,y,y’,y’’,…….,y(n-1)) (2a) với các điều kiện đầu: y(i)(x0) = yi0 i = 0,1,2,….,(n-1) (2b) Trong đó ta ký hiệu y(i) là đạo hàm cấp i của y hàm số y = y(x) Dễ dàng chuyển hệ (1) sang hệ 2 bằng cách đưa ra: y1= y, y2= y’, y3= y’’,…..,yn= y(n-1) Do đó ta có hệ phương trình gồm n phương trình vi phân cấp 1: y1’ = y2 y2’ = y3 …….. yn-1’ = yn yn’ = f(x,y1,y2,…..,yn) Với các điều kiện đầu: y1(x0) = y0, y2(x0) = y10,……,yn(x0)=yn-1,0 Sau đó ta áp dụng phương pháp RUNGE_KUTTA đối với hệ phương trình vi phân cấp 1 để giải phương trình trên. II/ Phương trình vi phân chuyển động của con lắc eliptic và cách giải bằng phương pháp RUNGE_KUTTA Con lắc eliptic là một cơ hệ gồm một lò xo có độ cứng c nối với vật có khối lượng m1. Vật có khối lượng m1 được nối với một chất điểm có khối lượng m2 không đáng kể thông qua thanh dẫn có chiều dài l. Động năng của con lắc có dạng: T = T1 + T2 Trong đó: T1: Động năng của vật có khối lượng m1 T2: Động năng của con lắc khối lượng m2 T1 = m1v12 = m12 T2 = m2v22= m2(12 + 12 ) x1 = lcos(j) ị 1 = -lsin(j) y1 = y + lsin(j) ị 1 = + lcos(j) ị (12 + 12 ) = l2sin2(j)2 + 2 + 2lcos(j) + l2cos2(j)2 = l22 + 2 + 2lcos(j) ị T= (m1 + m2) 2 + m2 (l22 + 2lcos(j)) Thế năng của cơ hệ: P = -m2glcos(j) + cy2 Theo định lý LAGRANGE II ta có: - = i= 1,…,n = (m1 + m2) + m2 lcos(j) = (m1 + m2) - m2 lsin(j)2 + m2 lcos(j) = 0,= cy = m2l2+ m2 lcos(j) = m2l2- m2lsin(j)+m2 lcos(j) = - m2lsin(j) = m2glsin(j) Hệ phương trình vi phân chuyển động của con lắc có dạng: (m1+m2) + m2lcos(j) - m2l2sin(j) + cy = 0 m2lcos(j) + m2l2+ m2glsin(j) = 0 Sử dụng các phép biến đổi toán học đưa hệ phương trình trên về dạng = (m2gsin(j)cos(j) + m2lsin(j)2- cx)/(m1 + m2sin2(j)) = (-cos(j)- gsin(j))/l Để thuận tiện cho việc lập trình đặt ff[1]=(m2gsin(y[3])cos(y[3])+m2lsin(y[3])y[4]y[4]–cy[1])/(m1+m2sin2(y[3]) ) ff[2] = (-cos(y[3])ff[1]-gsin(y[3])/l f[1] = y[2] f[2] = ff[1] f[3] = y[4] f[4] = ff[2] III/ Phần chương trình Phần này nhằm giải hệ phương trình vi phân cấp một và cấp 2 qua đó khảo sát cơ hệ con lắc eliptic . Chương trình được viết trên ngôn ngữ C. #include "fstream.h" #include"stdio.h" #include"math.h" #include"conio.h" #include"graphics.h" #include"dos.h" #include"iostream.h //****khoi tao do hoa****// void khoitao() { int mh=0,mode=0; initgraph(&mh,&mode,"d:\\tc\\bgi"); } //****Module goi ham thu nhat****// void dif(float x,float*y,float*f) { f[1]=0*x+y[2]; f[2]=(1-y[1]*y[1])*y[2]-y[1]; } //*****Module goi ham thu hai*****// void module2(float*ff,float*f,float*y,float x,float m1,float m2,float l,float g,float c) { ff[1]=(m2*g*cos(y[3])*sin(y[3])+m2*l*sin(y[3])*y[4]*y[4]-c*y[1])/(m1+m2*sin(y[3])*sin(y[3])); ff[2]=(0*x-g*sin(y[3])-cos(y[3])*ff[1])/l; f[1]=y[2]; f[2]=ff[1]; f[3]=y[4]; f[4]=ff[2]; } //*****Ham lam tron so*****// int round(float a) { if ((a-floor(a))>= 0.5) a=ceil(a); if ((a-floor(a))< 0.5) a=floor(a); return a; } //****Module giai he phuong trinh vi phan cap 1****// void runge(float x,float h,int nstep,int n,float*y,float*f, float*f1,float*f2,float*f3,float*f4,float*yy) { float ts,hh,tt,t1; float p3[200],p4[200][4]; int i,j,xd,yd; hh=0.5*h; ts=x; printf("bien x cac gia tri cua y"); for(i=1;i<=nstep;++i) { printf("\n %8.2f",x); p3[i]=x; for(j=1;j<=n;++j) printf("%16.2f",y[j]); for(j=1;j<=n;++j) p4[i][j]=y[j]; dif(x,y,f); for(j=1;j<=n;++j) f1[j]=h*f[j]; tt=x+hh; for(j=1;j<=n;++j) yy[j]=y[j]+0.5*f1[j]; dif(tt,yy,f); for(j=1;j<=n;++j) { f2[j]=h*f[j]; yy[j]=y[j]+0.5*f2[j]; } dif(tt,yy,f); tt=x+h; for(j=1;j<=n;++j) { f3[j]=h*f[j]; yy[j]=y[j]+f3[j]; } dif(tt,yy,f); x=ts+h*i; for(j=1;j<=n;++j) { f4[j]=h*f[j]; y[j]=y[j]+(f1[j]+2*f2[j]+2*f3[j]+f4[j])/6; } if (i%22==0) { getch();clrscr(); } } printf("\n Nhan ENTER de xem do thi cua co he"); getch(); //***Doan chuong trinh ve do thi cua co he***// khoitao(); setbkcolor(15); setcolor(1); rectangle(4,4,633,475); setcolor(8); xd=getmaxx()/2;yd=getmaxy()/2; line(xd-290,yd,xd+290,yd); line(xd,yd+200,xd,yd-180); line(xd+283,yd+3,xd+290,yd); line(xd+283,yd-3,xd+290,yd); line(xd+3,yd-173,xd,yd-180); line(xd-3,yd-173,xd,yd-180); setcolor(5); outtextxy(xd+150,yd+150,"DRAWING.."); setcolor(1); outtextxy(xd-8,yd+5,"O"); outtextxy(xd-30,yd-180,"x,v"); outtextxy(xd+288,yd+5,"t"); for(j=1;j<=n;++j) for(i=1;i<=nstep;++i) { if(j==n) setlinestyle(1,1,1); else setlinestyle(0,1,1); if(i==1) moveto(xd+round(25*p3[i]),yd-round(25*p4[i][j])); else lineto(xd+round(25*p3[i]),yd-round(25*p4[i][j])); delay(20); } setcolor(1); setlinestyle(0,1,1); line(xd+200,yd-155,xd+220,yd-155); setlinestyle(3,1,1); line(xd+200,yd-140,xd+220,yd-140); setcolor(0); outtextxy(xd+150,yd+150,"DRAWING.."); setcolor(1); outtextxy(xd+225,yd-160,"Vitri"); outtextxy(xd+225,yd-145,"Vantoc"); outtextxy(510,10,"PHAM TUAN LONG"); outtextxy(510,25,"CO TIN A-K43"); outtextxy(xd-120,yd-200,"DO THI MO TA DAO DONG CUA CO HE"); getch(); closegraph(); } //****Module giai bai toan ap dung cho he co hoc****// void runge2(float x,float h,int nstep,int m,float*y,float*f,float*ff, float*f1,float*f2,float*f3,float*f4,float*yy, float m1,float m2,float l,float g,float c) { float ts,hh,tt,xd,yd; float p1[200],p2[200][4]; int i,j,ii; char*dso[]={"0","10","20","30","40","50","60","70","80","90","100","110"}; char*dso1[]={"0","1","2","3","4","5"}; char*dso2[]={"0","-1","-2","-3","-4","-5"}; char*dso3[]={"0","-10","-20","-30","-40","-50","-60","-70","-80"}; hh=0.5*h; ts=x; ofstream kq("Ql"); printf(" bien x cac gia tri cua y"); for(i=1;i<=nstep;++i) { printf("\n %8.2f",x); p1[i]=x; kq<<x<<" "; for(j=1;j<=m;++j) { if(j==3) { y[j]=y[j]*M_PI/180; } printf("%12.2f",y[j]); kq<<y[j]<< " "; } kq<<endl; for(j=1;j<=m;++j)p2[i][j]=y[j]; module2(ff,f,y,x,m1,m2,l,g,c); for(j=1;j<=m;++j) f1[j]=h*f[j]; tt=x+hh; for(j=1;j<=m;++j) yy[j]=y[j]+0.5*f1[j]; module2(ff,f,yy,tt,m1,m2,l,g,c); for(j=1;j<=m;++j) { f2[j]=h*f[j]; yy[j]=y[j]+0.5*f2[j]; } module2(ff,f,yy,tt,m1,m2,l,g,c); tt=x+h; for(j=1;j<=m;++j) { f3[j]=h*f[j]; yy[j]=y[j]+f3[j]; } module2(ff,f,yy,tt,m1,m2,l,g,c); for(j=1;j<=m;++j) { f4[j]=h*f[j]; if(j==3) y[j]=(y[j]+(f1[j]+2*f2[j]+2*f3[j]+f4[j])/6)*180/M_PI; else y[j]=y[j]+(f1[j]+2*f2[j]+2*f3[j]+f4[j])/6; } x=ts+h*i; if (i%22==0) { getch();clrscr(); } } printf("\n Nhan ENTER de xem do thi cua co he"); getch(); //***Doan chuong trinh ve do thi cua co he***// khoitao(); setbkcolor(15); setcolor(1); rectangle(4,4,633,475); setcolor(8); xd=getmaxx()/2-100;yd=getmaxy()/2; line(xd-200,yd,xd+340,yd); line(xd,yd+180,xd,yd-180); line(xd+333,yd+3,xd+340,yd); line(xd+333,yd-3,xd+340,yd); line(xd+3,yd-173,xd,yd-180); line(xd-3,yd-173,xd,yd-180); setcolor(1); for(ii=1;ii<=11;ii++) { line(xd+30*ii,yd-2,xd+30*ii,yd+2); outtextxy(xd+30*ii-5,yd+4,dso[ii]); } for(ii=1;ii<=5;ii++) { line(xd+2,yd-30*ii,xd-2,yd-30*ii); outtextxy(xd-9,yd-30*ii-2,dso1[ii]); } for(ii=1;ii<=5;ii++) line(xd+2,yd+30*ii,xd-2,yd+30*ii); for(ii=1;ii<=5;ii++) outtextxy(xd-18,yd+30*ii-4,dso2[ii]); for(ii=1;ii<=6;ii++) line(xd-30*ii,yd-2,xd-30*ii,yd+2); for(ii=1;ii<=8;ii++) outtextxy(xd-30*ii-14,yd+4,dso3[ii]); setcolor(5); outtextxy(xd+200,yd+150,"DRAWING.."); setcolor(1); outtextxy(xd-8,yd+5,"O"); outtextxy(xd-30,yd-180,"x,v"); outtextxy(xd+340,yd-8,"t(s)"); setlinestyle(0,1,1); for(j=1;j<=m;++j) for(i=1;i<=nstep;++i) { if(j==2) setlinestyle(1,1,1); if(j==3) { setcolor(4); setlinestyle(0,1,1); } if(j==4) { setcolor(4); setlinestyle(1,1,1); } if(i==1) moveto(xd+round(30*p1[i]),yd-round(30*p2[i][j])); else lineto(xd+round(30*p1[i]),yd-round(30*p2[i][j])); delay(20); } setcolor(1); setlinestyle(0,1,1); line(xd+300,yd-185,xd+320,yd-185); setlinestyle(2,1,1); line(xd+300,yd-170,xd+320,yd-170); setcolor(4); setlinestyle(0,1,1); line(xd+300,yd-155,xd+320,yd-155); setlinestyle(2,1,1); line(xd+300,yd-140,xd+320,yd-140); setcolor(0); outtextxy(xd+200,yd+150,"DRAWING.."); setcolor(1); outtextxy(xd+330,yd-190,"y(t)"); outtextxy(xd+330,yd-175,"dy(t)"); setcolor(4); outtextxy(xd+330,yd-160,"phi(t)"); outtextxy(xd+330,yd-145,"dphi(t)"); setcolor(1); outtextxy(510,10,"PHAM TUAN LONG"); outtextxy(510,25,"CO TIN A-K43"); outtextxy(xd-120,yd-200,"DO THI MO TA DAO DONG CUA CO HE"); getch(); closegraph(); } //********chuong trinh chinh*********// void main() { float x,h,ch,m1,m2,l,g,c; float y[10],yy[10],f1[10],f2[10],f3[10],f4[10],f[10],ff[10]; int ma,n,nstep,m,i; tt:clrscr(); printf(" ******************CHUONG_TRINH_GIAI_PTVT_THEO_RUNGEKUTTA********************\n"); printf(" ============================================================================"); n=2;m=4;g=9.8; printf("\n\nChon1:Giai he ptvp VANDERPOL theo RUNGE_KUTTA"); printf("\nChon2:Ptvp cap 2 cho bai toan con lac ELIPTIC "); printf("\nBan muon giai he ptvp cap 1 hay 2:");scanf("%d",&ma); switch(ma) { case 1: printf("Nhap so doan can tinh nstep:");scanf("%d",&nstep); printf("Nhap buoc can tinh h:");scanf("%f",&h); printf("Cho gia tri ban dau cua x:");scanf("%f",&x); for(i=1;i<=n;++i) { printf("Cho gia tri ban dau cua y[%d]:",i);scanf("%f",&y[i]); } runge(x,h,nstep,n,y,f,f1,f2,f3,f4,yy); break; case 2: printf("Nhap so doan can tinh nstep:");scanf("%d",&nstep); printf("Nhap buoc can tinh h:");scanf("%f",&h); printf("Cho gia tri ban dau cua x:");scanf("%f",&x); printf("Cho khoi luong cua vat 1 m1:");scanf("%f",&m1); printf("Cho khoi luong cua vat 2 m2:");scanf("%f",&m2); printf("Cho chieu dai l:");scanf("%f",&l); printf("Cho do cung cua lo xo c:");scanf("%f",&c); for(i=1;i<=m;++i) { printf("Cho gia tri ban dau cua y[%d]:",i);scanf("%f",&y[i]); } runge2(x,h,nstep,m,y,f,ff,f1,f2,f3,f4,yy,m1,m2,l,g,c); break; } printf("\n Ban co muon tinh tiep tuc nua khong(C/K)?"); ch=getch(); if(ch=='c'||ch=='C') goto tt; } Trong đó: nstep: Số các bước cần tính h: Khoảng cách giữa 2 nút kế tiếp m1: Khối lượng của vật một m2: Khối lượng của con lắc l: Chiều dài thanh dẫn c: Độ cứng của lò xo y[1]: Vị trí ban đầu của vật có khối lượng m1 y[2]: Vận tốc ban đầu của vật có khối lượng m1 y[3]: Vị trí ban đầu của con lắc y[4]: Vận tốc ban đầu của con lắc Đó là chương trình giải hệ phương trình vi phân cấp một và hai bằng phương pháp RUNGE_KUTTA. Molule 1 nhằm giải hệ phương trình vi phân cấp một, module hai nhằm giải hệ phương trình vi phân cấp hai áp dụng cho bài toán cơ. Khi cần sử dụng chương trình với các hệ phương trình vi phân khác nhau người sử dụng chỉ cần thay đổi các module sau: Đối với hệ phương trình vi phân cấp một chỉ cần tác động vào module: void dif(float x,float*y,float*f) { …………………… …………………… } Người sử dụng chỉ việc đưa số phương trình cần tính toán vào module trên và thay đổi số n ở chương trình chính (n cho biết số phương trình cần đưa và để tính toán ) Đối với hệ phương trình vi phân cấp hai, người sử dụng chỉ cần tác động vào module: void module2(float*ff,float*f,float*y,float x,float m1,float m2,float l,float g,float c) { ………………………………………… ………………………………………… } Sau đó thay đổi số phương trình cần tính toán tại chương trình chính. void main() { ……………………………………… ……………………………………… } IV/ Một số kết quả của chương trình Với nstep =110 ; m1= 50; m2= 4; l=10; c=100 y[1] = 2; y[2]=0; y[3] = 0; y[4] = 0; Ta có đồ thị: Với nstep =110 ; m1= 50; m2= 4; l=10; c=100 y[1] = 0; y[2]=2; y[3] = 0; y[4] = 0; Ta có đồ thị mô tả dao động của con lắc. Với nstep =110 ; m1=50; m2=4; l=10; c=100 y[1] = 0; y[2]=0; y[3]=60o; y[4]=0; Ta có đồ thị: Với nstep =110 ; m1=50; m2=4; l=10; c=100 y[1] = 0; y[2]=0; y[3]=0; y[4]=2; Ta có đồ thị: Kết luận: Với phương pháp RUNGE_KUTTA, ta có thể giải được các phương trình phi tuyến với độ chính xác khá cao và có thể khảo sát chuyển động của những cơ hệ phức tạp qua đó nghiên cứu được đặc tính của cơ hệ như vị trí,vận tốc, gia tốc… Tài liệu tham khảo Đinh Văn Phong: Phương pháp số trong cơ học. NXB Khoa Học Kỹ Thuật. Tạ Văn Đĩnh: Phương pháp tính. NXB Giáo Dục. Đỗ Sanh: Cơ lý thuyết (Tĩnh học và động học). NXB Giáo Dục. Đỗ Sanh: Cơ lý thuyết (Động lực học). NXB Giáo Dục. Mục Lục Trang Lời nói đầu 1 Phần thuyết minh 2 I/ Giải gần đúng hệ phương trình vi phân thường bằng 2 Phương pháp RUNGE_KUTTA II/ Phương trình vi phân chuyển động của con lắc 3 eliptic và cách giải bằng phương pháp RUNGE_KUTTA: III/ Phần chương trình 5 IV/ Một số kết quả của chương trình 15 ._.