Numerical method in python- can't spot the problem?

I am writing this numerical method formula of trapezium rule for double integrals. enter image description here

Note that hx = (b-a)/nx, hy = (d-c)/ny to get the interval widths and xj = a+hxj and yi = c+hyi

edited Nov 23 '18 at 15:39

asked Nov 22 '18 at 22:46

NewYork

243

1

Why did you remove the code?

– kvantour
Nov 23 '18 at 16:32

I didn't think it made sense

– NewYork
Mar 4 at 18:42

add a comment |

I am writing this numerical method formula of trapezium rule for double integrals. enter image description here

Note that hx = (b-a)/nx, hy = (d-c)/ny to get the interval widths and xj = a+hxj and yi = c+hyi

edited Nov 23 '18 at 15:39

asked Nov 22 '18 at 22:46

NewYork

243

1

Why did you remove the code?

– kvantour
Nov 23 '18 at 16:32

I didn't think it made sense

– NewYork
Mar 4 at 18:42

add a comment |

I am writing this numerical method formula of trapezium rule for double integrals. enter image description here

Note that hx = (b-a)/nx, hy = (d-c)/ny to get the interval widths and xj = a+hxj and yi = c+hyi

edited Nov 23 '18 at 15:39

asked Nov 22 '18 at 22:46

NewYork

243

I am writing this numerical method formula of trapezium rule for double integrals. enter image description here

Note that hx = (b-a)/nx, hy = (d-c)/ny to get the interval widths and xj = a+hxj and yi = c+hyi

python python-3.x math integration numerical-methods

edited Nov 23 '18 at 15:39

asked Nov 22 '18 at 22:46

NewYork

243

edited Nov 23 '18 at 15:39

asked Nov 22 '18 at 22:46

NewYork

243

edited Nov 23 '18 at 15:39

asked Nov 22 '18 at 22:46

NewYork

243

asked Nov 22 '18 at 22:46

NewYork

243

asked Nov 22 '18 at 22:46

NewYork

243

1

Why did you remove the code?

– kvantour
Nov 23 '18 at 16:32

I didn't think it made sense

– NewYork
Mar 4 at 18:42

add a comment |

1

Why did you remove the code?

– kvantour
Nov 23 '18 at 16:32

I didn't think it made sense

– NewYork
Mar 4 at 18:42

Why did you remove the code?

– kvantour
Nov 23 '18 at 16:32

I didn't think it made sense

– NewYork
Mar 4 at 18:42

add a comment |

1 Answer
1

active

oldest

votes

A few problems in your code:

First yes your indentation here is off (but I assume it's from not copying it across well since this would lead to an error rather than a wrong value). In the future make sure the indentation in your question corresponds to what you have at on your own computer before posting...

Then a term should be added within a for if and only if it's in the corresponding sum... Here you put everything within the double for loop which corresponds to having all the terms in the double sum.

Finally range(1,n) already stops at n-1 only so you want to remove those -1 in the ranges.

In the end:

def double_integral(f,a,b,c,d,nx,ny):



    hx = (b-a)/nx

    hy = (d-c)/ny 



    first_term = (f(a,c)+f(a,d)+f(b,c)+f(b,d))



    i_sum = 0

    for i in range(1,ny):

        i_sum += f(a,c+i*hy)+f(b, c+i*hy)



    j_sum = 0

    for j in range(1,nx):

        j_sum += f(a+j*hx,c)+f(a+j*hx,d)



    ij_sum = 0

    for i in range(1,ny):

        for j in range(1,nx):

            ij_sum += f(a+j*hx,c+i*hy)



    integral = (first_term/4 + i_sum/2 + j_sum/2 + ij_sum) * hx * hy



    return integral





def t(x,y):

    return x*(y**(2))



print(double_integral(t,0,2,0,1,10,10))



0.6700000000000003

You'll get closer to 2/3 by choosing numbers of steps larger than 10...

And you can be more pythonic by using sum comprehension:

def double_integral(f,a,b,c,d,nx,ny):

    hx = (b-a)/nx

    hy = (d-c)/ny 

    first_term = (f(a,c)+f(a,d)+f(b,c)+f(b,d))

    i_sum = sum(f(a,c+i*hy)+f(b, c+i*hy) for i in range (1,ny))

    j_sum = sum(f(a+j*hx,c)+f(a+j*hx,d) for j in range(1,nx))

    ij_sum = sum(f(a+j*hx,c+i*hy) for i in range (1,ny) for j in range(1,nx))

    integral = (first_term/4 + i_sum/2 + j_sum/2 + ij_sum) * hx * hy

    return integral

edited Nov 22 '18 at 23:18

answered Nov 22 '18 at 23:08

Julien

7,70831637

Thank you so much. This is really helpful!

– NewYork
Nov 23 '18 at 1:14

Can I ask what do you mean by your second point: 'Then a term should be added within a for if and only if it's in the corresponding sum... Here you put everything within the double for loop which corresponds to having all the terms in the double sum.' Is it invalid to write for i in range:... and then for j in range....?

– NewYork
Nov 23 '18 at 1:21

1

If a term is inside both for i and for j loops it will be added for all combinations of values of i and j: sum(i for i in range(3)) == 0+1+2 == 3 but sum(i for j in range(10) for i in range(3)) == 30 because you are adding each value of i over and over 10 times for each value of j.

– Julien
Nov 23 '18 at 2:32

add a comment |

Your Answer

StackExchange.ifUsing("editor", function () {
StackExchange.using("externalEditor", function () {
StackExchange.using("snippets", function () {
StackExchange.snippets.init();
});
});
}, "code-snippets");

StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "1"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});

function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: true,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});

}
});

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53438773%2fnumerical-method-in-python-cant-spot-the-problem%23new-answer', 'question_page');
}
);

Post as a guest

Name

Required, but never shown

1 Answer
1

active

oldest

votes

1 Answer
1

active

oldest

votes

A few problems in your code:

Finally range(1,n) already stops at n-1 only so you want to remove those -1 in the ranges.

In the end:

def double_integral(f,a,b,c,d,nx,ny):



    hx = (b-a)/nx

    hy = (d-c)/ny 



    first_term = (f(a,c)+f(a,d)+f(b,c)+f(b,d))



    i_sum = 0

    for i in range(1,ny):

        i_sum += f(a,c+i*hy)+f(b, c+i*hy)



    j_sum = 0

    for j in range(1,nx):

        j_sum += f(a+j*hx,c)+f(a+j*hx,d)



    ij_sum = 0

    for i in range(1,ny):

        for j in range(1,nx):

            ij_sum += f(a+j*hx,c+i*hy)



    integral = (first_term/4 + i_sum/2 + j_sum/2 + ij_sum) * hx * hy



    return integral





def t(x,y):

    return x*(y**(2))



print(double_integral(t,0,2,0,1,10,10))



0.6700000000000003

You'll get closer to 2/3 by choosing numbers of steps larger than 10...

And you can be more pythonic by using sum comprehension:

def double_integral(f,a,b,c,d,nx,ny):

    hx = (b-a)/nx

    hy = (d-c)/ny 

    first_term = (f(a,c)+f(a,d)+f(b,c)+f(b,d))

    i_sum = sum(f(a,c+i*hy)+f(b, c+i*hy) for i in range (1,ny))

    j_sum = sum(f(a+j*hx,c)+f(a+j*hx,d) for j in range(1,nx))

    ij_sum = sum(f(a+j*hx,c+i*hy) for i in range (1,ny) for j in range(1,nx))

    integral = (first_term/4 + i_sum/2 + j_sum/2 + ij_sum) * hx * hy

    return integral

edited Nov 22 '18 at 23:18

answered Nov 22 '18 at 23:08

Julien

7,70831637

Thank you so much. This is really helpful!

– NewYork
Nov 23 '18 at 1:14

Can I ask what do you mean by your second point: 'Then a term should be added within a for if and only if it's in the corresponding sum... Here you put everything within the double for loop which corresponds to having all the terms in the double sum.' Is it invalid to write for i in range:... and then for j in range....?

– NewYork
Nov 23 '18 at 1:21

1

If a term is inside both for i and for j loops it will be added for all combinations of values of i and j: sum(i for i in range(3)) == 0+1+2 == 3 but sum(i for j in range(10) for i in range(3)) == 30 because you are adding each value of i over and over 10 times for each value of j.

– Julien
Nov 23 '18 at 2:32

add a comment |

A few problems in your code:

Finally range(1,n) already stops at n-1 only so you want to remove those -1 in the ranges.

In the end:

def double_integral(f,a,b,c,d,nx,ny):



    hx = (b-a)/nx

    hy = (d-c)/ny 



    first_term = (f(a,c)+f(a,d)+f(b,c)+f(b,d))



    i_sum = 0

    for i in range(1,ny):

        i_sum += f(a,c+i*hy)+f(b, c+i*hy)



    j_sum = 0

    for j in range(1,nx):

        j_sum += f(a+j*hx,c)+f(a+j*hx,d)



    ij_sum = 0

    for i in range(1,ny):

        for j in range(1,nx):

            ij_sum += f(a+j*hx,c+i*hy)



    integral = (first_term/4 + i_sum/2 + j_sum/2 + ij_sum) * hx * hy



    return integral





def t(x,y):

    return x*(y**(2))



print(double_integral(t,0,2,0,1,10,10))



0.6700000000000003

You'll get closer to 2/3 by choosing numbers of steps larger than 10...

And you can be more pythonic by using sum comprehension:

def double_integral(f,a,b,c,d,nx,ny):

    hx = (b-a)/nx

    hy = (d-c)/ny 

    first_term = (f(a,c)+f(a,d)+f(b,c)+f(b,d))

    i_sum = sum(f(a,c+i*hy)+f(b, c+i*hy) for i in range (1,ny))

    j_sum = sum(f(a+j*hx,c)+f(a+j*hx,d) for j in range(1,nx))

    ij_sum = sum(f(a+j*hx,c+i*hy) for i in range (1,ny) for j in range(1,nx))

    integral = (first_term/4 + i_sum/2 + j_sum/2 + ij_sum) * hx * hy

    return integral

edited Nov 22 '18 at 23:18

answered Nov 22 '18 at 23:08

Julien

7,70831637

Thank you so much. This is really helpful!

– NewYork
Nov 23 '18 at 1:14

Can I ask what do you mean by your second point: 'Then a term should be added within a for if and only if it's in the corresponding sum... Here you put everything within the double for loop which corresponds to having all the terms in the double sum.' Is it invalid to write for i in range:... and then for j in range....?

– NewYork
Nov 23 '18 at 1:21

1

If a term is inside both for i and for j loops it will be added for all combinations of values of i and j: sum(i for i in range(3)) == 0+1+2 == 3 but sum(i for j in range(10) for i in range(3)) == 30 because you are adding each value of i over and over 10 times for each value of j.

– Julien
Nov 23 '18 at 2:32

add a comment |

A few problems in your code:

Finally range(1,n) already stops at n-1 only so you want to remove those -1 in the ranges.

In the end:

def double_integral(f,a,b,c,d,nx,ny):



    hx = (b-a)/nx

    hy = (d-c)/ny 



    first_term = (f(a,c)+f(a,d)+f(b,c)+f(b,d))



    i_sum = 0

    for i in range(1,ny):

        i_sum += f(a,c+i*hy)+f(b, c+i*hy)



    j_sum = 0

    for j in range(1,nx):

        j_sum += f(a+j*hx,c)+f(a+j*hx,d)



    ij_sum = 0

    for i in range(1,ny):

        for j in range(1,nx):

            ij_sum += f(a+j*hx,c+i*hy)



    integral = (first_term/4 + i_sum/2 + j_sum/2 + ij_sum) * hx * hy



    return integral





def t(x,y):

    return x*(y**(2))



print(double_integral(t,0,2,0,1,10,10))



0.6700000000000003

You'll get closer to 2/3 by choosing numbers of steps larger than 10...

And you can be more pythonic by using sum comprehension:

def double_integral(f,a,b,c,d,nx,ny):

    hx = (b-a)/nx

    hy = (d-c)/ny 

    first_term = (f(a,c)+f(a,d)+f(b,c)+f(b,d))

    i_sum = sum(f(a,c+i*hy)+f(b, c+i*hy) for i in range (1,ny))

    j_sum = sum(f(a+j*hx,c)+f(a+j*hx,d) for j in range(1,nx))

    ij_sum = sum(f(a+j*hx,c+i*hy) for i in range (1,ny) for j in range(1,nx))

    integral = (first_term/4 + i_sum/2 + j_sum/2 + ij_sum) * hx * hy

    return integral

edited Nov 22 '18 at 23:18

answered Nov 22 '18 at 23:08

Julien

7,70831637

A few problems in your code:

Finally range(1,n) already stops at n-1 only so you want to remove those -1 in the ranges.

In the end:

def double_integral(f,a,b,c,d,nx,ny):



    hx = (b-a)/nx

    hy = (d-c)/ny 



    first_term = (f(a,c)+f(a,d)+f(b,c)+f(b,d))



    i_sum = 0

    for i in range(1,ny):

        i_sum += f(a,c+i*hy)+f(b, c+i*hy)



    j_sum = 0

    for j in range(1,nx):

        j_sum += f(a+j*hx,c)+f(a+j*hx,d)



    ij_sum = 0

    for i in range(1,ny):

        for j in range(1,nx):

            ij_sum += f(a+j*hx,c+i*hy)



    integral = (first_term/4 + i_sum/2 + j_sum/2 + ij_sum) * hx * hy



    return integral





def t(x,y):

    return x*(y**(2))



print(double_integral(t,0,2,0,1,10,10))



0.6700000000000003

You'll get closer to 2/3 by choosing numbers of steps larger than 10...

And you can be more pythonic by using sum comprehension:

def double_integral(f,a,b,c,d,nx,ny):

    hx = (b-a)/nx

    hy = (d-c)/ny 

    first_term = (f(a,c)+f(a,d)+f(b,c)+f(b,d))

    i_sum = sum(f(a,c+i*hy)+f(b, c+i*hy) for i in range (1,ny))

    j_sum = sum(f(a+j*hx,c)+f(a+j*hx,d) for j in range(1,nx))

    ij_sum = sum(f(a+j*hx,c+i*hy) for i in range (1,ny) for j in range(1,nx))

    integral = (first_term/4 + i_sum/2 + j_sum/2 + ij_sum) * hx * hy

    return integral

edited Nov 22 '18 at 23:18

answered Nov 22 '18 at 23:08

Julien

7,70831637

edited Nov 22 '18 at 23:18

answered Nov 22 '18 at 23:08

Julien

7,70831637

answered Nov 22 '18 at 23:08

Julien

7,70831637

answered Nov 22 '18 at 23:08

Julien

7,70831637

Thank you so much. This is really helpful!

– NewYork
Nov 23 '18 at 1:14

Can I ask what do you mean by your second point: 'Then a term should be added within a for if and only if it's in the corresponding sum... Here you put everything within the double for loop which corresponds to having all the terms in the double sum.' Is it invalid to write for i in range:... and then for j in range....?

– NewYork
Nov 23 '18 at 1:21

1

If a term is inside both for i and for j loops it will be added for all combinations of values of i and j: sum(i for i in range(3)) == 0+1+2 == 3 but sum(i for j in range(10) for i in range(3)) == 30 because you are adding each value of i over and over 10 times for each value of j.

– Julien
Nov 23 '18 at 2:32

add a comment |

Thank you so much. This is really helpful!

– NewYork
Nov 23 '18 at 1:14

Can I ask what do you mean by your second point: 'Then a term should be added within a for if and only if it's in the corresponding sum... Here you put everything within the double for loop which corresponds to having all the terms in the double sum.' Is it invalid to write for i in range:... and then for j in range....?

– NewYork
Nov 23 '18 at 1:21

1

If a term is inside both for i and for j loops it will be added for all combinations of values of i and j: sum(i for i in range(3)) == 0+1+2 == 3 but sum(i for j in range(10) for i in range(3)) == 30 because you are adding each value of i over and over 10 times for each value of j.

– Julien
Nov 23 '18 at 2:32

Thank you so much. This is really helpful!

– NewYork
Nov 23 '18 at 1:14

Can I ask what do you mean by your second point: 'Then a term should be added within a for if and only if it's in the corresponding sum... Here you put everything within the double for loop which corresponds to having all the terms in the double sum.' Is it invalid to write for i in range:... and then for j in range....?

– NewYork
Nov 23 '18 at 1:21

If a term is inside both for i and for j loops it will be added for all combinations of values of i and j: sum(i for i in range(3)) == 0+1+2 == 3 but sum(i for j in range(10) for i in range(3)) == 30 because you are adding each value of i over and over 10 times for each value of j.

– Julien
Nov 23 '18 at 2:32

add a comment |

draft saved

draft discarded

Thanks for contributing an answer to Stack Overflow!

Please be sure to answer the question. Provide details and share your research!

But avoid …

Asking for help, clarification, or responding to other answers.

Making statements based on opinion; back them up with references or personal experience.

To learn more, see our tips on writing great answers.

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

Post as a guest

Name

Required, but never shown

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

Name

Required, but never shown

Name

Required, but never shown

This page is only for reference, If you need detailed information, please check here

搜尋此網誌

Argthtjtr