Kết quả tính tích phân bằng Monte-Carlo không chính xác

#1

Mọi người xem giúp em bài này code đã đúng chưa với ạ. Đề bài là: “Tính tích phân từ 1 đến 2 của hàm f(x) =x^2”. Em tính bằng máy thì ra 2,3333333 nhưng mà code xong thì kết quả toàn < 2 thôi ạ.

Em cảm ơn mọi người ạ !!

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#define MAX 99

float f(float x)
{
    return (x*x);
}
int main()
{
    int N, m=0, n=0,i=0;
    double S;
    float x[99999];
    float y[99999];
	printf("Nhap so diem muon tao: ");
	scanf("%d", &N);
	for( ; i<N; i++)
    {
        x[i]=rand() %2 + 1;
        y[i]=rand() %4 + 1;
        if(y[i]>f(x[i]))
            {
                m++;
            }
        else
        {
            n++;
        }
    }
    printf("%d", n);
    S = (double)((2-1)*(4-1)*(n))/(N);
    printf("\nTich phan tu 1 den 2 cua f(x) = x^2 la: %10.6f", S);
	return 0;
}
#2

Dùng rand() thì thêm câu lệnh srand(time(NULL); vào đầu chương trình nhé bạn :slight_smile:

#3

e thêm rồi nhưng nó vẫn ra kết quả sai bác ạ !!

#4

Mình cũng chưa hiểu lắm về cách làm của bạn. Nhưng nếu dùng toán thì hàm ban đầu chỉ là 1/3 * x^3 còn nếu dùng rand() kia thì kết quả sẽ có nhiều biến động và làm sao cố định được?

#5

rand() toàn ra số nguyên =))

Đây gọi là cách tính diện tích kiểu phóng phi tiêu (thuật ngữ gọi là Monte-Carlo :p), để tính tích phân gần đúng có thể dùng hình thang để xấp xỉ cái hình (xem tài liệu môn pp tính).

Chỉ cần lưu seed thôi.

4 Likes
#6

yêu cầu chỉ cần tính gần đúng thôi ạ. tức là kiểu cho càng nhiều điểm thì càng chính xác đó ạ !! nhưng vấn đề em thấy kết quả nó ngáo ngáo

#7

hàm ban đầu mà là 1/3*x^3 thì nó là nguyên hàm của x^2 luôn rồi bác :)) em đang muốn làm từ hàm đơn giản ý. vì có thể đôi hàm không biểu diễn được nguyên hàm mà

#8

Lỗi chỗ này

S = (double)((2-1)*(4-1)*(n))/(N);

Bạn thay (4-1) chỉ còn 4 thôi, được:

S = (double)((2-1)*4*(n))/(N);

Vì đoạn code trừ đi 1 nên nó mất đi phần diện tích ứng với 1 <= x <= 2, 0 <= y <= 1. Trong khi diện tích đó thuộc tích phân.

3 Likes
#9

mình thử thay rồi nhưng kết quả vẫn sai bạn ạ

#10

Đương nhiên nó sẽ có sai số lớn, mình chỉ chỉnh lại công thức cho chính xác.
Còn để ra kết quả chính xác tới 7/3 thì có bạn @rog10 đã gợi ý cho bạn.

"rand() trả về số nguyên", nghĩa là x thuộc { 1; 2 } và y thuộc { 1; 2; 3; 4 }, nên (x;y) có 8 cặp giá trị: { (1, 1); (1, 2), (1, 3), (1, 4), (2, 1), (2, 2), (2, 3), (2, 4) }.

Các giá trị thuộc diện tích có: (1, 1), (2, 1), (2, 2), (2, 3), (2, 4), nên tỉ số các điểm trong tích phân và tổng số điểm là 5/8

Diện tích giới hạn 1 <= x <= 2, 0 <= y <= 4 là S = 1*4 = 4

Tích phân tính xấp xỉ: 5/8 * 4 = 2.5
Do đó khi chạy chương trình chỉ xấp xỉ quanh 2.5

Chỉ có 8 giá trị (x;y) đó khó mà ra kết quả chính xác được. :sweat_smile:

4 Likes
#11

Đừng dùng float, sai số lớn lắm, dùng double ấy

#12

em cảm ơn bác ạ, em hiểu rồi ạ

#13

cho ra số nguyên từ 1, 2 hoặc 1, 2, 3, 4 chứ ko phải số thực :V

với lại tích phân của x^2 là dưới x^2 nên hình vuông random phải là từ (1,0) tới (2,4) chứ random trong (1,1) tới (2,4) thì thiếu mất phần (1,0) tới (2,1) có diện tích đúng bằng 1 nên lúc nào cũng < 2 là đúng rồi :V nếu +1 vào thì nó ra > 2 thôi

3 Likes
#14 
        x[i]= (double) (rand()%2 + 1);
        y[i]= (double) (rand()%5 + 0);

Như này thì nó ra số thực đúng ko ạ?

#15

ko :V là giá trị nguyên kiểu số thực: 0.0, 1.0, 2.0, 3.0, 4.0

2 Likes
#16

thế muốn tạo ra số thực thì làm sao ạ?
em mới học C nên hơi mơ hồ ạ

#17

Mình thì mình sẽ ko dùng rand vì nó rất ẹ về mọi nhẽ, nhưng đại loại là ntn: rand() cho max từ 0 đến RAND_MAX, vậy làm sao cho nó từ 0 đến 1 :smiley:

2 Likes
#18 
double randf() // số thực ngẫu nhiên trong [0, 1]
{
    return rand() % RAND_MAX;
}

double randf_range(double a, double b) // số thực ngẫu nhiên trong [a, b]
{
    return randf() * (b - a) + a;
}

x[i] = randf_range(1, 2);
y[i] = randf_range(0, 4);
3 Likes
#19

C++ hiện đại :V

#include <chrono>
#include <iostream>
#include <random>
using namespace std::chrono;

double f(double x)
{
    return x * x;
}

int main()
{
    // nếu xài MSVC 2017 thì lấy auto seed = std::random_device{}();
    // nếu có xài Boost thì lấy auto seed = boost::random::random_device{}();
    // nếu muốn đơn giản hơn thì cứ lấy auto seed = time(0); nhớ thêm <ctime>
    auto seed =
        duration_cast<seconds>(system_clock::now().time_since_epoch()).count();
    std::mt19937 prbg{seed};
    int cnt = 0;
    int ptCnt = 100000;
    for (int i = 0; i < ptCnt; ++i)
    {
        auto x = std::uniform_real_distribution<>{1.0, 2.0}(prbg);
        auto y = std::uniform_real_distribution<>{0.0, 4.0}(prbg);
        if (y <= f(x)) cnt++;
    }
    std::cout << (2.0 - 1.0) * (4.0 - 0.0) * cnt / ptCnt << "\n";
}
4 Likes
83% thành viên diễn đàn không hỏi bài tập, còn bạn thì sao?