..::Vật Lý Việt Nam::..  


Trang Chủ Tạp Chí VLVN Thư Viện Thành Viên Xuất Sắc Bảng Xếp Hạng Theo Tháng Bản Rút Gọn
Trở lại   ..::Vật Lý Việt Nam::.. > Modern Physics- Vật lý hiện đại > Computational Physics and Computer Science
Ghi Danh Hỏi/ÐápLuật Diễn Đàn Thành Viên Lịch Tìm Kiếm Bài Trong Ngày Ðánh Dấu Ðã Ðọc

Trả lời
 
Ðiều Chỉnh Xếp Bài
Old 26-09-2010, 11:49 PM   #1
hunter800
Junior Member
 
Tham gia: Sep 2010
Quốc gia:
Giới tính: Male
Bài gửi: 2
hunter800 is on a distinguished road
Chê bai: 0
Bị chê 0 lần trong 0 bài
Cám ơn: 0
Được cám ơn 0 lần trong 0 bài
Unhappy Giải hệ phương trình vi phân nhiều biến trong Matlab

Mình đang làm 1 project để mô phỏng trường điện từ truyền trong ống dẫn sóng.Bây h cần xây dựng 6 phương trình vi phân 3 biến của 6 thành phần trường.Cho mình ỏi Matlab có cho xây dựng và giải hệ phương trình vi phân nhiều biến không?Xin các bạn giúp mình với!

thay đổi nội dung bởi: pththao, 27-09-2010 lúc 12:20 AM. Lý do: Sửa chính tả tên đề tài.
hunter800 is offline   Trả Lời Với Trích Dẫn
Old 27-09-2010, 12:08 AM   #2
pththao
forget me not
 
pththao's Avatar
 
Tham gia: Apr 2010
Cư trú: casablanca
Quốc gia:
Giới tính: Male
Bài gửi: 183
pththao is on a distinguished road
Chê bai: 2
Bị chê 0 lần trong 0 bài
Cám ơn: 95
Được cám ơn 135 lần trong 85 bài
Theo hiểu biết của mình xây dựng phương trình thì nên làm bằng tay. Hoặc nhác có thể dùng symbolic để thiết lập, nhưng người viết phải khá khéo léo để chủ động thời gian và chuyển đổi các dạng dữ liệu. Việc giải phương trình vi phân nhiều biến là "nói chung" có thể thực hiện, tất nhiên là chủ yếu numerically: Việc này có thể làm bằng lệnh ODE23 hoặc ODE45.
pththao is offline   Trả Lời Với Trích Dẫn
Old 30-09-2010, 11:34 PM   #3
hunter800
Junior Member
 
Tham gia: Sep 2010
Quốc gia:
Giới tính: Male
Bài gửi: 2
hunter800 is on a distinguished road
Chê bai: 0
Bị chê 0 lần trong 0 bài
Cám ơn: 0
Được cám ơn 0 lần trong 0 bài
Unhappy

bạn có thể nói rõ hơn được không?
ví dụ giải t vi phân sau trong matlab thì làm ntn?

y=a*dy/dx+b*dy/dt+c*dy/dz

y1=e*dy1/dz+ f*dy/dz+ g*dy1/dt

.....

thank bạn nhiều
hunter800 is offline   Trả Lời Với Trích Dẫn
Old 24-11-2010, 12:00 AM   #4
anh012anh
Junior Member
 
Tham gia: Oct 2008
Quốc gia:
Giới tính: Male
Bài gửi: 1
anh012anh is on a distinguished road
Chê bai: 0
Bị chê 0 lần trong 0 bài
Cám ơn: 0
Được cám ơn 0 lần trong 0 bài
Giai he phuong trinh vi phan cap 1

Hien jo minh dang lam luan van cao hoc,su dung matlab,minh dang "bi" khi viet chuong trinh de giai he phuong trinh vi phan cap 1.
Vi du:
x ' = a11 (t ) x + a12 (t ) y + f1 (t )

⎩ y ' = a21 (t ) x + a22 (t ) y + f 2 (t )
Cac ban nao gioi matlab cho i kien voi....thanks
anh012anh is offline   Trả Lời Với Trích Dẫn
Old 24-12-2010, 11:21 PM   #5
vllt_Hue
Administrator
 
Tham gia: Dec 2010
Quốc gia:
Giới tính: Male
Bài gửi: 362
vllt_Hue is on a distinguished road
Chê bai: 0
Bị chê 2 lần trong 2 bài
Cám ơn: 1
Được cám ơn 454 lần trong 255 bài
Giái số hệ phương trình vi phân trong Matlab

Để giải số các phương trình và hệ phương trình vi phân, thuật toán thông dụng nhât là Runge _ Kutta bậc 4 (RK4) (Nếu muốn chi tiết bạn phải xem các sách về Phương pháp số). Ơ đây minh đưa ra đoạn chương trình trên Matlab để giải hệ hai phương trình vi phân câp 1 theo thuat toan tren :
x' = a11(t) x + a12(t) y + f1(t)
y' = a21(t) x + a12(t) y + f2(t)
với dieu kien dau x (t0) = x0; y(t0) = y0 và t0<= t <= t1 (Các phương pháp giải số phai cho biết điều kiện đầu và khoảng cần tính)
%---- Day la dong chu thich trong Matlab
clear all;
% Cac dieu kien dau
x0 = ;% nhap các gia trị của bài toán cụ thể vào day
y0 = ;
t0 = ;
t1 = ;
dt = ; % đây là bước để tính số, càng nhỏ kết quả sẽ chính xác hơn
t = t0:dt:t1;
n = length(t);
x = zeros(n,1); % các mảng chứa kết quả tính số
y = zeros(n,1);
x(1) = x0;
y(1) = y0;
dx = @(x,y,t) a11(t)*x+a12(t)*y+f1(t) ;% Nhập biểu thức cụ thể vào đây
dy = @(x,y,t) a21(t)*x+a22(t)*y+f2(t) ;
% Cách khai bao hàm kiểu này dùng cho Matlab 7 trở lên
for i = 1:n-1
k1x = dt*dx(x(i),y(i),t(i));
k1y = dt*dy(x(i),y(i),t(i));
k2x = dt*dx(x(i)+k1x/2,y(i)+k1y/2,t(i)+dt/2);
k2y = dt*dy(x(i)+k1x/2,y(i)+k1y/2,t(i)+dt/2);
k3x = dt*dx(x(i)+k2x/2,y(i)+k2y/2,t(i)+dt/2);
k3y = dt*dy(x(i)+k2x/2,y(i)+k2y/2,t(i)+dt/2);
k4x = dt*dx(x(i)+k3x,y(i)+k3y,t(i)+dt);
k4y = dt*dy(x(i)+k3x,y(i)+k3y,t(i)+dt);
x(i+1) = x(i) + (k1x + 2*k2x + 2*k3x + k4x)/6;
y(i+1) = y(i) + (k1y + 2*k2y + 2*k3y + k4y)/6;
end
% Vẽ đồ thị kết quả tính số
% Vẽ đồ thị x theo t
subplot(4,1,1)
plot(t,x)
% Vẽ đồ thị y theo t
subplot(4,1,2)
plot(t,y)
% Vẽ đồ thi của x, y theo t trên cùng 1 trục
subplot(4,1,3)
plot(t,x,t,y)
% Ve đồ thi của y theo x
subplot(4,1,4)
plot(x,y)
% Neu muon lưu kết qua vào file data thì làm như sau
% Tao mang gồm n hang 3 cot, cot 1 chứa t, cột 2 chưa x, cột 3 chứa y
data = zeros(n,3)
data(:,1) = t;
data(:,2) = x;
data(:,3) = y;
% Đặt tên file lưu dữ liệu
fna = 'C:\ketqua.dat';
% Lưu dữ liệu từ mang vào file
save(fna,'data','-ASCII')

% Khi có file dữ liệu rồi bạn có thể sử dụng phần mềm Origin để vẽ đồ thị
% Phân mềm nay ở trong này cũng có giới thiệu rồi đó

thay đổi nội dung bởi: vllt_Hue, 24-12-2010 lúc 11:31 PM.
vllt_Hue is offline   Trả Lời Với Trích Dẫn
Các thành viên gửi lời cám ơn tới vllt_Hue vì bài viết này:
day.dreamer (25-12-2010), quockhanh185 (02-08-2011), Vạn lý Độc hành (25-12-2010)
Old 25-12-2010, 10:41 AM   #6
nhatbmt771
Junior Member
 
Tham gia: Dec 2010
Quốc gia:
Giới tính: Male
Bài gửi: 7
nhatbmt771 is on a distinguished road
Chê bai: 0
Bị chê 0 lần trong 0 bài
Cám ơn: 4
Được cám ơn 2 lần trong 2 bài
@vllt_Hue: Còn giải hệ phương trình vi phân bằng giải thuật của Euler thì sao hả bác:
Ví dụ như hệ này:


Giải thuật của Euler thế này:
nhatbmt771 is offline   Trả Lời Với Trích Dẫn
Các thành viên gửi lời cám ơn tới nhatbmt771 vì bài viết này:
Vạn lý Độc hành (25-12-2010)
Old 26-12-2010, 09:52 PM   #7
vllt_Hue
Administrator
 
Tham gia: Dec 2010
Quốc gia:
Giới tính: Male
Bài gửi: 362
vllt_Hue is on a distinguished road
Chê bai: 0
Bị chê 2 lần trong 2 bài
Cám ơn: 1
Được cám ơn 454 lần trong 255 bài
Cách giải gần giống như trên (thuật toán Euler có độ chính xác kém hơn RK4 nhưng đơn giản hơn) với mỗi hàm y ta cần một hàm dy để biểu diễn vế phải của phương trình và hai biến trung gian K1,K2 .
Chuong trinh chi tiêt như sau :
clear all;
% Cac dieu kien dau
x0 = ;
x1 = ;
y10 = ;
y20 = ;
y30 = ;
% Buoc tinh
dx = ;
x = x0:dx:x1;
n = length(x);
y1 = zeros(n,1);
y2 = zeros(n,1);
y3 = zeros(n,1);
y1(1) = y10;
y2(1) = y20;
y3(1) = y30;
dy1 = @(x,u1,u2,u3) u2 ;
dy2 = @(x,u1,u2,u3) u3 ;
dy3 = @(x,u1,u2,u3) 2*u3+2*u2 -2*u1 +(sin(x)+4*x)*exp(2*x);
for i = 1:n-1
% Tinh k1 cho y1, y2, y3
k11 = dx*dy1(x(i),y1(i),y2(i),y3(i));
k12 = dx*dy2(x(i),y1(i),y2(i),y3(i));
k13 = dx*dy3(x(i),y1(i),y2(i),y3(i));
% Tinh k2 cho y1, y2, y3
k21 = dt*dy1(x(i) + dx,y1(i) + k11, y2(i) + k12, y3(i) + k13);
k22 = dt*dy2(x(i) + dx,y1(i) + k11, y2(i) + k12, y3(i) + k13);
k23 = dt*dy3(x(i) + dx,y1(i) + k11, y2(i) + k12, y3(i) + k13);
% cap nhat gia tri y
y1(i+1) = y1(i) + (k11 + k21)/2;
y2(i+1) = y2(i) + (k12 + k22)/2;
y3(i+1) = y2(i) + (k13 + k23)/2;
end
% Ở đây ta giả một phương trình vi phân cấp 3 bằng cách đưa về hệ 3 phương trình vi phân cấp 1, do đó kết quả là y3
% Ve đồ thị
plot(x,y3)
% luu ket qua
data = zeros(n,2)
data(:,1) = x;
data(:,2) = y3;
% Đặt tên file lưu dữ liệu
fna = 'C:\ketqua.dat';
% Lưu dữ liệu từ mang vào file
save(fna,'data','-ASCII')
vllt_Hue is offline   Trả Lời Với Trích Dẫn
Các thành viên gửi lời cám ơn tới vllt_Hue vì bài viết này:
nhatbmt771 (27-12-2010), Vạn lý Độc hành (27-12-2010)
Old 27-12-2010, 03:58 PM   #8
nhatbmt771
Junior Member
 
Tham gia: Dec 2010
Quốc gia:
Giới tính: Male
Bài gửi: 7
nhatbmt771 is on a distinguished road
Chê bai: 0
Bị chê 0 lần trong 0 bài
Cám ơn: 4
Được cám ơn 2 lần trong 2 bài
@vllt_HUE: Cám ơn bạn đã trả lời.
Mình có vài thắc mắc sau:
Thứ nhất là đề mình đưa bị nhầm một xíu là ở phương trình y3' sửa 2y2 thành y2. Mà cái này thì ko quan trọng lắm.

Trên code của bạn còn thiếu định nghĩa cho dt. Và theo euler thì nó bằng dx luôn, mình đã bổ sung.
Trích:
clear all;
% Cac dieu kien dau
x0 =0 ;
x1 =2 ;
y10 =1 ;
y20 =2 ;
y30 =3 ;
% Buoc tinh
dx =0.1 ;dt =0.1 ;
x = x0:dx:x1;
n = length(x);
y1 = zeros(n,1);
y2 = zeros(n,1);
y3 = zeros(n,1);
y1(1) = y10;
y2(1) = y20;
y3(1) = y30;
dy1 = @(x,u1,u2,u3) u2 ;
dy2 = @(x,u1,u2,u3) u3 ;
dy3 = @(x,u1,u2,u3) 2*u3+u2 -2*u1 +(sin(x)+4*x)*exp(2*x);
for i = 1:n-1
% Tinh k1 cho y1, y2, y3
k11 = dx*dy1(x(i),y1(i),y2(i),y3(i));
k12 = dx*dy2(x(i),y1(i),y2(i),y3(i));
k13 = dx*dy3(x(i),y1(i),y2(i),y3(i));
% Tinh k2 cho y1, y2, y3
k21 = dt*dy1(x(i) + dx,y1(i) + k11, y2(i) + k12, y3(i) + k13);
k22 = dt*dy2(x(i) + dx,y1(i) + k11, y2(i) + k12, y3(i) + k13);
k23 = dt*dy3(x(i) + dx,y1(i) + k11, y2(i) + k12, y3(i) + k13);
% cap nhat gia tri y
y1(i+1) = y1(i) + (k11 + k21)/2;
y2(i+1) = y2(i) + (k12 + k22)/2;
y3(i+1) = y2(i) + (k13 + k23)/2;
end
% Ở đây ta giả một phương trình vi phân cấp 3 bằng cách đưa về hệ 3 phương trình vi phân cấp 1, do đó kết quả là y3
% Ve đồ thị
plot(x,y3)
% luu ket qua
data = zeros(n,2)
data(:,1) = x;
data(:,2) = y3;
% Đặt tên file lưu dữ liệu
fna = 'C:\ketqua.dat';
% Lưu dữ liệu từ mang vào file
save(fna,'data','-ASCII')

Và đây là kết quả sau khi mình sửa thành y2 và thêm dt, được save tại fide ketqua.dat
Trích:
0.0000000e+000 3.0000000e+000
1.0000000e-001 2.6855249e+000
2.0000000e-001 3.0110650e+000
3.0000000e-001 3.4723729e+000
4.0000000e-001 4.0451874e+000
5.0000000e-001 4.7551019e+000
6.0000000e-001 5.6401014e+000
7.0000000e-001 6.7492677e+000
8.0000000e-001 8.1453203e+000
9.0000000e-001 9.9081833e+000
1.0000000e+000 1.2139490e+001
1.1000000e+000 1.4968227e+001
1.2000000e+000 1.8557793e+001
1.3000000e+000 2.3114818e+001
1.4000000e+000 2.8900184e+001
1.5000000e+000 3.6242771e+001
1.6000000e+000 4.5556609e+001
1.7000000e+000 5.7362271e+001
1.8000000e+000 7.2313537e+001
1.9000000e+000 9.1230641e+001
2.0000000e+000 1.1514170e+002

có phải 1.1514170e+002 là 115.14170 không ạ?

Và mấy ngày nay mình cũng tự mày mò, thử đủ kiểu thì mình tự code hàm mình như sau:

Trích:
function [t,u,w,x]=Euler(f,g,z,x0,t0,u0,w0,h,x1)
if nargin<4, error('Ham co it nhat 4 doi so');end;
x=[];x=[x,x0];n=(x1-x0)/h;
for i=1:n, x(i+1)=x(i)+h;end;
t=[];t=[t,t0];
u=[];u=[u,u0];
w=[];w=[w,w0];
for i=1:n
K1A=h*f(x(i),t(i),u(i),w(i));
K1B=h*g(x(i),t(i),u(i),w(i));
K1C=h*z(x(i),t(i),u(i),w(i));
K2A=h*f(x(i)+h,t(i)+K1A,u(i)+K1B,w(i)+K1C);
K2B=h*g(x(i)+h,t(i)+K1A,u(i)+K1B,w(i)+K1C);
K2C=h*z(x(i)+h,t(i)+K1A,u(i)+K1B,w(i)+K1C);
t(i+1)=t(i)+0.5*(K1A+K2A);
u(i+1)=u(i)+0.5*(K1B+K2B);
w(i+1)=w(i)+0.5*(K1C+K2C);
end;

sau đó mình nhập các giá trị ban đầu, nhập hàm thông qua lệnh inline.
Trích:
>> x0=0;x1=2;h=0.1;t0=1;u0=2;w0=3;
>> f=inline('u','x','t','u','w')
f =
Inline function:
f(x,t,u,w) = u
>> g=inline('w','x','t','u','w');
>> z=inline('2*w+u-2*t+(sin(x)+4*x)*exp(2*x)','x','t','u','w');

Và cho ra kết quả là.
Trích:
t =
Columns 1 through 9
1.0000 1.2150 1.4664 1.7632 2.1171 2.5443 3.0662 3.7120 4.5215
Columns 10 through 18
5.5481 6.8645 8.5686 10.7921 13.7121 17.5652 22.6675 29.4395 38.4388
Columns 19 through 21
50.4032 66.3057 87.4260
u =
Columns 1 through 9
2.0000 2.3300 2.7380 3.2495 3.8999 4.7376 5.8285 7.2626 9.1617
Columns 10 through 18
11.6901 15.0686 19.5929 25.6572 33.7857 44.6723 59.2335 78.6755 104.5817
Columns 19 through 21
139.0247 184.7109 245.1659
w =
Columns 1 through 9
3.0000 3.6855 4.5916 5.8015 7.4288 9.6271 12.6038 16.6369 22.0979
Columns 10 through 18
29.4814 39.4429 52.8483 70.8379 94.9082 127.0183 169.7254 226.3596 301.2467
Columns 19 through 21
399.9944 529.8573 700.2024
x =
Columns 1 through 9
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000
Columns 10 through 18
0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000
Columns 19 through 21
1.8000 1.9000 2.0000

So kết quả của mình là w tương ứng với y3 bên code của bạn. Nó khác nhau quá :(.

Vậy nhờ bạn đánh giá giùm mình cái code, với kết quả trên. Cám ơn bạn rất nhiều.

thay đổi nội dung bởi: nhatbmt771, 27-12-2010 lúc 04:11 PM.
nhatbmt771 is offline   Trả Lời Với Trích Dẫn
Old 27-12-2010, 04:24 PM   #9
nhatbmt771
Junior Member
 
Tham gia: Dec 2010
Quốc gia:
Giới tính: Male
Bài gửi: 7
nhatbmt771 is on a distinguished road
Chê bai: 0
Bị chê 0 lần trong 0 bài
Cám ơn: 4
Được cám ơn 2 lần trong 2 bài
Trời ơi.! Đã tìm ra chỗ sai rồi.
Bạn code sai ở một chỗ chết người là:
Trích:
y3(i+1) = y2(i) + (k13 + k23)/2;

Đúng ra là y3 chứ không phải là y2 :eek:

Giờ thì kết quả là đúng rồi. Cám ơn bạn rất nhiều, mình biết thêm 1 cách code nhờ bạn .

ah, với lại cho mình hỏi là giờ nếu muốn lấy kết quả sau dấu , nhiều hơn 4 chữ số thì chỉnh ở đâu ạ?
nhatbmt771 is offline   Trả Lời Với Trích Dẫn
Các thành viên gửi lời cám ơn tới nhatbmt771 vì bài viết này:
Vạn lý Độc hành (27-12-2010)
Old 27-12-2010, 10:24 PM   #10
vllt_Hue
Administrator
 
Tham gia: Dec 2010
Quốc gia:
Giới tính: Male
Bài gửi: 362
vllt_Hue is on a distinguished road
Chê bai: 0
Bị chê 2 lần trong 2 bài
Cám ơn: 1
Được cám ơn 454 lần trong 255 bài
+ Mình bị nhầm ở chỗ các công thức tính k21,k22,k23 vì ở đây giải phương trình của y theo x, do đó bước tính là dx chứ không phải dt. Do đó trong đoạn code của mình bạn thay dt trong đó bởi dx là được;
+ 1.1514170e+002 = 1.1514170*10^2= 115.14170
+ Trong Matlab có hai định dạng khi hiển thị các số thực là shortlong. Ngầm định là short với 4 chữ số lẻ sau dấu phẩy. Để thay đổi, ở đầu chương trình bạn dùng lệnh format long là được. Khi đó bạn sẽ được 14 đến 15 chữ số lẻ sau dấu phẩy
vllt_Hue is offline   Trả Lời Với Trích Dẫn
Các thành viên gửi lời cám ơn tới vllt_Hue vì bài viết này:
nhatbmt771 (28-12-2010), Vạn lý Độc hành (27-12-2010)
Old 27-12-2010, 11:22 PM   #11
nhatbmt771
Junior Member
 
Tham gia: Dec 2010
Quốc gia:
Giới tính: Male
Bài gửi: 7
nhatbmt771 is on a distinguished road
Chê bai: 0
Bị chê 0 lần trong 0 bài
Cám ơn: 4
Được cám ơn 2 lần trong 2 bài
Vậy phương trình của mình là đúng chứ?
Với lại nghiệm cuối là w vậy thì theo pt gốc nó tương đương với y'' chứ đâu phải y''' cần tìm?

Với lại cho mình hỏi đề có 1 câu là: Mô tả bằng đồ thị so với nghiệm chính xác của phương trình.

Vậy thì phải hiểu và làm như thế nào?
nhatbmt771 is offline   Trả Lời Với Trích Dẫn
Old 28-12-2010, 01:57 AM   #12
vllt_Hue
Administrator
 
Tham gia: Dec 2010
Quốc gia:
Giới tính: Male
Bài gửi: 362
vllt_Hue is on a distinguished road
Chê bai: 0
Bị chê 2 lần trong 2 bài
Cám ơn: 1
Được cám ơn 454 lần trong 255 bài
Post

Code cua ban không sai, tuy nhiên có một số chỗ hơi thừa
function [t,u,w,x]=Euler(f,g,z,x0,t0,u0,w0,h,x1)
if nargin<4, error('Ham co it nhat 4 doi so');end;
% Đoạn này bạn chỉ cần khai bao nhu sau
x = x0:h:x1 ; % nhu vay là duoc mảng các giá trị của x từ x0 -> x1, bước h
n = length(x) : % tra về số phần tử của mang x;
% Không cần phải dùng vòng lặp for để tính x(i)
x=[];x=[x,x0];n=(x1-x0)/h;
for i=1:n, x(i+1)=x(i)+h;end;

t=[];t=[t,t0];
u=[];u=[u,u0];
w=[];w=[w,w0];
for i=1:n
K1A=h*f(x(i),t(i),u(i),w(i));
K1B=h*g(x(i),t(i),u(i),w(i));
K1C=h*z(x(i),t(i),u(i),w(i));
K2A=h*f(x(i)+h,t(i)+K1A,u(i)+K1B,w(i)+K1C);
K2B=h*g(x(i)+h,t(i)+K1A,u(i)+K1B,w(i)+K1C);
K2C=h*z(x(i)+h,t(i)+K1A,u(i)+K1B,w(i)+K1C);
t(i+1)=t(i)+0.5*(K1A+K2A);
u(i+1)=u(i)+0.5*(K1B+K2B);
w(i+1)=w(i)+0.5*(K1C+K2C);
end;
* Mình bị nhầm, với hệ phương trình của bạn thì hàm cần tìm là y1(x) = y(x) tức là t trong hàm của bạn
* Mô tả bằng đồ thị so với nghiệm chính xác của phương trình :
Ở đây yêu cầu so sánh kết quả giải số và kết quả giải chính xác bằng cách vẽ đồ thị nghiệm gần đúng và nghiệm chính xác trên cùng một trục.
Phương trình gốc của bạn là :
y'''(x)-2y''(x)-y'(x)+2y(x)=(sin(x) + 4x)exp(2x)
Phương trình này có thể giải chính xác bằng phương pháp hệ số bất định của lý thuyết phương trình vi phân. trong Matlab bạn có thể giải chính xác nó bằng lệnh dsolve như sau : Nhập dòng lệnh sau vào cửa sổ lệnh của Matlab :
dsolve('D3y-2*D2y-Dy+2*y=(sin(x)+4*x)*exp(2*x)','x')
Trong nghiệm sẽ chứa 3 hằng số tùy ý được tìm từ điều kiện đầu.
Nếu có điều kiện đầu bạn có thể dùng lệnh :
dsolve('D3y-2*D2y-Dy+2*y=(sin(x)+4*x)*exp(2*x)','D2y(x)=.. ','Dy(x0)= ..','y(x0) = ','x')
Ở đây ban đưa các giá trị cụ thê của y(x0),y'(x0) và y''(x0) vào
vllt_Hue is offline   Trả Lời Với Trích Dẫn
Các thành viên gửi lời cám ơn tới vllt_Hue vì bài viết này:
Mai Thuận (01-05-2014), nhatbmt771 (28-12-2010)
Old 28-12-2010, 12:11 PM   #13
nhatbmt771
Junior Member
 
Tham gia: Dec 2010
Quốc gia:
Giới tính: Male
Bài gửi: 7
nhatbmt771 is on a distinguished road
Chê bai: 0
Bị chê 0 lần trong 0 bài
Cám ơn: 4
Được cám ơn 2 lần trong 2 bài
Trích:
Nguyên văn bởi vllt_Hue
* Mô tả bằng đồ thị so với nghiệm chính xác của phương trình :
Ở đây yêu cầu so sánh kết quả giải số và kết quả giải chính xác bằng cách vẽ đồ thị nghiệm gần đúng và nghiệm chính xác trên cùng một trục.
Phương trình gốc của bạn là :
y'''(x)-2y''(x)-y'(x)+2y(x)=(sin(x) + 4x)exp(2x)
Phương trình này có thể giải chính xác bằng phương pháp hệ số bất định của lý thuyết phương trình vi phân. trong Matlab bạn có thể giải chính xác nó bằng lệnh dsolve như sau : Nhập dòng lệnh sau vào cửa sổ lệnh của Matlab :
dsolve('D3y-2*D2y-Dy+2*y=(sin(x)+4*x)*exp(2*x)','x')
Trong nghiệm sẽ chứa 3 hằng số tùy ý được tìm từ điều kiện đầu.
Nếu có điều kiện đầu bạn có thể dùng lệnh :
dsolve('D3y-2*D2y-Dy+2*y=(sin(x)+4*x)*exp(2*x)','D2y(x)=.. ','Dy(x0)= ..','y(x0) = ','x')
Ở đây ban đưa các giá trị cụ thê của y(x0),y'(x0) và y''(x0) vào

Mình nhập hàm:
dsolve('D3y-2*D2y-Dy+2*y=(sin(x)+4*x)*exp(2*x)','D2y(x0)=3 ','Dy(x0)= 2','y(x0) = 1 ','x')

Và đây là kết quả được trả về:
Trích:
ans =
(exp(2*x)*(120*x - 9*cos(x) + 27*sin(x) - 40))/540 - exp(2*x)*(cos(x)/3 - (2*x^2)/3) - (exp(2*x)*(8*x - cos(x) + sin(x) - 8))/4 + (exp(2*x)*(exp(2*x0)*cos(x0) - 2*x0^2*exp(2*x0) + 2))/(3*exp(2*x0)) + (exp(x)*(8*x0*exp(2*x0) - 8*exp(2*x0) - exp(2*x0)*cos(x0) + exp(2*x0)*sin(x0) + 2))/(4*exp(x0)) - (exp(x0)*(120*x0*exp(2*x0) - 40*exp(2*x0) - 9*exp(2*x0)*cos(x0) + 27*exp(2*x0)*sin(x0) + 90))/(540*exp(x))

Vậy làm sao để vẽ hàm này lên đồ thị?
Mình click vào plot bên khung Workspace, chọn plot 2D. thì nó thực hiện hàm như và báo lỗi như sau:
Trích:
>> plot (ans, 'DisplayName', 'ans', 'YDataSource', 'ans'); figure(gcf)
??? Error using ==> plot
Conversion to double from sym is not possible.

Thử dùng lệnh như sau:
Trích:
>> y=dsolve('D3y-2*D2y-Dy+2*y=(sin(x)+4*x)*exp(2*x)','D2y(x0)=3 ','Dy(x0)= 2','y(x0) = 1 ','x');

>> x=linspace(0,2);
>> plot (x,y)
??? Error using ==> plot
Conversion to double from sym is not possible.

Bạn có thể nói thêm cụ thể về phần này không? trình tự vẽ 2 hàm trên cùng 1 đồ thị là như thế nào?
Cảm ơn bạn nhiều .

thay đổi nội dung bởi: nhatbmt771, 28-12-2010 lúc 12:31 PM.
nhatbmt771 is offline   Trả Lời Với Trích Dẫn
Old 28-12-2010, 01:28 PM   #14
vllt_Hue
Administrator
 
Tham gia: Dec 2010
Quốc gia:
Giới tính: Male
Bài gửi: 362
vllt_Hue is on a distinguished road
Chê bai: 0
Bị chê 2 lần trong 2 bài
Cám ơn: 1
Được cám ơn 454 lần trong 255 bài
*Bạn phải thay các giá trị ban đầu của x0 vào nữa. Giả sử x0 = 0, khi đó các lệnh để giải và vẽ đồ thị là :
y = dsolve('D3y-2*D2y-Dy+2*y=(sin(x)+4*x)*exp(2*x)','D2y(0)=3','Dy(0)=2' ,'y(0)=1','x');
f = inline(y);
x=linspace(0,2);plot(x,f(x))
* Để so sánh với kết quả giải số, bạn giải số với cùng điều kiện đầu và khoảng ( ở đây là (0,2)) và giả sử kết quả tính số là t và u(t) Để so sánh đồ thị của hai kết quả bạn dùng lệnh :
plot(x,f(x),t,u)
vllt_Hue is offline   Trả Lời Với Trích Dẫn
Các thành viên gửi lời cám ơn tới vllt_Hue vì bài viết này:
nhatbmt771 (28-12-2010)
Old 28-12-2010, 01:49 PM   #15
nhatbmt771
Junior Member
 
Tham gia: Dec 2010
Quốc gia:
Giới tính: Male
Bài gửi: 7
nhatbmt771 is on a distinguished road
Chê bai: 0
Bị chê 0 lần trong 0 bài
Cám ơn: 4
Được cám ơn 2 lần trong 2 bài
Trích:
Nguyên văn bởi vllt_Hue
*Bạn phải thay các giá trị ban đầu của x0 vào nữa. Giả sử x0 = 0, khi đó các lệnh để giải và vẽ đồ thị là :
y = dsolve('D3y-2*D2y-Dy+2*y=(sin(x)+4*x)*exp(2*x)','D2y(0)=3','Dy(0)=2' ,'y(0)=1','x');
f = inline(y);
x=linspace(0,2);plot(x,f(x))
* Để so sánh với kết quả giải số, bạn giải số với cùng điều kiện đầu và khoảng ( ở đây là (0,2)) và giả sử kết quả tính số là t và u(t) Để so sánh đồ thị của hai kết quả bạn dùng lệnh :
plot(x,f(x),t,u)

Đoạn này mình không hiểu lắm.
Có phải là giả sử mình để nguyên hàm số mình tính thì đồ thị của pp số là (x,t)- hàm t theo x.
Rồi, đó là bước thứ nhất.
Bước 2 là hàm của bạn mình đổi lại là 1 hàm y theo z và đặt là hàm j (cho khỏi trùng với hàm f đã xài lúc tính bằng pp số), lúc đó mình cho cả 2 hàm số và hàm chính xác chạy. Sau đó: plot (z,j(z),x,t)
Mình thử mà vẫn báo lỗi...
nhatbmt771 is offline   Trả Lời Với Trích Dẫn
Old 28-12-2010, 01:54 PM   #16
nhatbmt771
Junior Member
 
Tham gia: Dec 2010
Quốc gia:
Giới tính: Male
Bài gửi: 7
nhatbmt771 is on a distinguished road
Chê bai: 0
Bị chê 0 lần trong 0 bài
Cám ơn: 4
Được cám ơn 2 lần trong 2 bài
àh, được rồi. Lý do là mình mới chỉ nhập hàm số chứ chưa cho chạy nó nên chương trình chưa load được giá trị miền t.
Vậy là ổn hết rồi, cám ơn bạn rất rất nhiều .
nhatbmt771 is offline   Trả Lời Với Trích Dẫn
Old 08-12-2011, 09:07 PM   #17
vungnguyengia
Junior Member
 
Tham gia: Dec 2011
Quốc gia:
Giới tính: Male
Bài gửi: 1
vungnguyengia is on a distinguished road
Chê bai: 0
Bị chê 0 lần trong 0 bài
Cám ơn: 0
Được cám ơn 0 lần trong 0 bài
x' = a11(t) x + a12(t) y + f1(t)
y' = a21(t) x + a12(t) y + f2(t)
Mọi người có thể cho mình đoạn code bài này giải bằng phương pháp khử không?
mình xin thanks trước.
vungnguyengia is offline   Trả Lời Với Trích Dẫn
Old 12-04-2013, 05:32 PM   #18
songlam0
Junior Member
 
Tham gia: Apr 2013
Quốc gia:
Giới tính: Male
Bài gửi: 3
songlam0 is on a distinguished road
Chê bai: 0
Bị chê 1 lần trong 1 bài
Cám ơn: 0
Được cám ơn 0 lần trong 0 bài
mìnhđang làm một đề tài với thầy giáo môn sức bền. khổ nỗi lòi ra hệ phuong trình vi phân cấp 2 mà ko gải bằng phương pháp bình thường được. lên mạng tìm nghe nói matlab có thể tìm ra được nghiệm gần đúng, thậm chí là chính xác. mún thử lắm. khổ nỗi thời gian thì gấp, cũng đang học matlab, nhưng chưa làm dc. các bác cao tay giúp em giải ra nghiệm với


A.d2y/dx+(M-n.x).dz/dx+(P-q.x).y=0
A.d2z/dx-(M-n.x).dy/dx+(P-q.x).z=0
tìm y(x) và z(x)
trong đó
(A,M,n,P,q là các hằng số
thực ra A là hạng tử tự do
M,P là mô men. n,q hình như là mô men gì đó, cũng ko rõ nữa)
các bác giải hộ em rùi cho em công thức nghiệm dc ko? em cần gấp lắm
songlam0 is offline   Trả Lời Với Trích Dẫn
Old 01-05-2014, 11:25 AM   #19
Mai Thuận
Junior Member
 
Tham gia: May 2014
Quốc gia:
Giới tính: Male
Bài gửi: 1
Mai Thuận is on a distinguished road
Chê bai: 0
Bị chê 0 lần trong 0 bài
Cám ơn: 1
Được cám ơn 0 lần trong 0 bài
Giải hệ phương trình vi phân bằng Mathlab

Hiện giờ mình đang làm luận văn cao học phần Quang, su dung matlab,minh dang bí khi viet chuong trinh de giai he phuong trinh vi phan:
Giả sử laser ngẫu nhiên phát hai mode có mật độ photon của các mode lần lượt là n1, n2. Hệ phương trình biểu diễn sự thay đổi mật độ photon của các mode theo thời gian có dạng:
dn1/dt = a1n1 - b1n1n1 - c12n1n2 + d12n2
dn2/dt = a2n2 - b2n2n2 - c21n2n1 + d21n1
Vẽ đồ thị biểu thị sự phụ thuộc của n1, n2 vào các tham số. Trong đó:
ai: là hệ số khuyếch đại.
bi: là hệ số mất mát.
cij: là hệ số liên kết trường.
dij: là hệ số photon hopping.

Sử dụng ngôn ngữ lập trình Matlab giải hệ phương trình tốc độ (2.1) và (2.2) để khảo sát hoạt động không dừng của laser ngẫu nhiên phát hai mode (Vẽ đồ thị biểu thị sự phụ thuộc của n1, n2 vào các tham số.)
Điều kiện ban đầu cho số photon của mode 1 và 2 tương ứng là n1(0) = 0.8 và n2(0) = 0.8 trong đơn vị tương đối (a.u).
Yêu cầu: Khảo sát xong thì đưa vào word cả 16 đồ thị và đồng thời gửi cả file Mathlab.
Copy cả bảng giá trị fu thuộc (n1, n2 fu thuộc vào 4 tham số).

thay đổi nội dung bởi: Mai Thuận, 01-05-2014 lúc 07:54 PM.
Mai Thuận is offline   Trả Lời Với Trích Dẫn
Old 09-05-2014, 10:02 PM   #20
vllt_Hue
Administrator
 
Tham gia: Dec 2010
Quốc gia:
Giới tính: Male
Bài gửi: 362
vllt_Hue is on a distinguished road
Chê bai: 0
Bị chê 2 lần trong 2 bài
Cám ơn: 1
Được cám ơn 454 lần trong 255 bài
Phương pháp giải ở đây giổng như bài toán mình đã trình bày trong bài viết ngày 24/12/2010, chỉ thay đổi một đôi chút trong hệ phương trình n1 thay cho x và n2 thay cho y. Mình viết lại cho bạn :
% đây là hệ phương trình của bạn
%dn1/dt = a1 n1 - b1*n1n1 - c12n1n2 + d12n2
%dn2/dt = a2n2 - b2n2n2 - c21n2n1 + d21n1
%Vẽ đồ thị biểu thị sự phụ thuộc của n1, n2 vào các tham số. Trong đó:
%ai: là hệ số khuyếch đại.
%bi: là hệ số mất mát.
%cij: là hệ số liên kết trường.
%dij: là hệ số photon hopping.
%---------------
clear all;
% Nhap cac tham so cua bai toan vao day
a1 = ;
b1 = ;
a2 = ;
b2 = ;
c12 = ;
c21 = ;
d12 = ;
d21 = ;
% Dieu kien dau (gia tri cua n1,n2 luc t = 0)
n10 = 0.8;
n20 = 0.8;
tmax = ; % Khoảng thời gian khảo sát
dt = ; % đây là bước để tính số, càng nhỏ kết quả sẽ chính xác hơn
t = 0:dt:tmax;
n = length(t);
x = zeros(n,1); % các mảng chứa kết quả tính số x-> n1, y-> n2
y = zeros(n,1);
x(1) = n10;
y(1) = n20;
dx = @(x,y,t) a1*x - b1*x^2 - c12*x*y + d12*y ;
dy = @(x,y,t) a2*y - b2*y^2 - c21*x*y + d21*x;
% Cách khai bao hàm kiểu này dùng cho Matlab 7 trở lên
for i = 1:n-1
k1x = dt*dx(x(i),y(i),t(i));
k1y = dt*dy(x(i),y(i),t(i));
k2x = dt*dx(x(i)+k1x/2,y(i)+k1y/2,t(i)+dt/2);
k2y = dt*dy(x(i)+k1x/2,y(i)+k1y/2,t(i)+dt/2);
k3x = dt*dx(x(i)+k2x/2,y(i)+k2y/2,t(i)+dt/2);
k3y = dt*dy(x(i)+k2x/2,y(i)+k2y/2,t(i)+dt/2);
k4x = dt*dx(x(i)+k3x,y(i)+k3y,t(i)+dt);
k4y = dt*dy(x(i)+k3x,y(i)+k3y,t(i)+dt);
x(i+1) = x(i) + (k1x + 2*k2x + 2*k3x + k4x)/6;
y(i+1) = y(i) + (k1y + 2*k2y + 2*k3y + k4y)/6;
end
% Vẽ đồ thị kết quả tính số
% Vẽ đồ thị x theo t
figure (1)
plot(t,x)
% Vẽ đồ thị y theo t
figure(2)
plot(t,y)
% Neu muon lưu kết qua vào file data thì làm như sau
% Tao mang gồm n hang 3 cot, cot 1 chứa t, cột 2 chưa x, cột 3 chứa y
data = zeros(n,3)
data(:,1) = t; %Cot 1 là thời gian
data(:,2) = x; % cot 2 là n1
data(:,3) = y; %cot 3 la n2
% Đặt tên file lưu dữ liệu
fna = 'C:\ketqua.dat';
% Lưu dữ liệu từ mang vào file
save(fna,'data','-ASCII')
%-----------------
Theo cách này với mỗi bộ giá trị của các tham số ta sẽ nhận được đồ thị phụ thuộc của n1 và n2 theo t.
vllt_Hue is offline   Trả Lời Với Trích Dẫn
Trả lời


Ðang đọc: 1 (0 thành viên và 1 khách)
 
Ðiều Chỉnh
Xếp Bài

Quyền Sử Dụng Ở Diễn Ðàn
Bạn không được quyền gởi bài
Bạn không được quyền gởi trả lời
Bạn không được quyền gởi kèm file
Bạn không được quyền sửa bài

vB code đang Mở
Smilies đang Mở
[IMG] đang Mở
HTML đang Tắt
Chuyển đến

Chủ đề tương tự
Ðề tài Người gửi Forum Trả lời Bài viết mới
Thuyết mới về vật lý học "Vật lý : Hạt -Trường thống nhất" duytho 1001 câu chuyện liên quan đến Vật lý 2 01-02-2010 12:26 PM
Nhìn lại những giải Nobel về Vật lí 1901-2004 you know who Lịch sử và danh nhân Vật lý 13 06-07-2006 12:11 AM




Powered by vBulletin, PhysicsVN Community
Copyright ©2014 All Rights Reserved.