Category Archives: Mathematica

math, Math Geek, Mathematica

A Barnsley’s Fern In 7 Lines of Mathematica

Leave a comment
I used to draw a Barnsley’s fern with a program written in Pascal when I was 19 years old. Yesterday someone asked about it in a forum I visited, so I drew another one using Mathematica. The code is much shorter this time. (I’m sure that many people can shorten it even more.)
Here’s the code to draw the fern with 10,000 points. You can copy and paste and run it in Mathematica:

          
ifsFern[p_] := Module[{i},
  i = Random[Integer, 99];
  If[i < 1, Return[{{0., 0.}, {0., .16}}.p ]];
  If[i >= 1 && i < 86, Return[{{0.85, 0.04}, {-0.04, 0.85}}.p + {0., 1.6}]];
  If[i >= 86 && i < 93, Return[{{0.20, -0.26}, {0.23, 0.22}}.p + {0., 1.6}]];
  If[i >= 93, Return[{{-0.15, 0.28}, {0.26, 0.24}}.p + {0., 0.44}]]] 

Graphics[{RGBColor[0, 0.5, 0], Point[NestList[ifsFern, {0, 0}, 10000]]}]
   





The result looks like this:
math, Mathematica, science, science class, ภาษาไทย, สอนเด็กๆ

แก้สมการด้วยวิธีของนิวตัน

1 Comment

 
มีเด็กๆที่สนใจคณิตศาสตร์มาถามผมว่าสมการทั่วๆไปแก้ออกมาเป๊ะๆไม่ได้แล้วเราทำอย่างไร ผมก็บอกว่าสมการส่วนใหญ่เราต้องหาคำตอบด้วยการประมาณเอาครับ ซึ่งวิธีอันหนึ่งที่เราสามารถใช้ได้ง่ายๆก็คือวิธีของนิวตัน ข้างล่างนี้เป็นกระทู้ที่ผมเขียนไว้ที่ Mahidol Physics Educational Center ครับ:

******

สำหรับปัญหาที่เราแก้สมการโดยตรงไม่ได้ เราต้องแก้ด้วยวิธีประมาณด้วยตัวเลขครับ วิธีที่ใช้กันบ่อยๆวิธีหนึ่งก็คือวิธีการของนิวตัน (Newton’s method: http://en.wikipedia.org/wiki/Newton’s_method) ครับ

วิธีการของนิวตันบอกว่า ถ้าจะแก้สมการ f(x) = 0 ให้เราเดาค่า x มาสักค่า (เรียกมันว่า x0) ก็แล้วกัน แล้วค่า x อันต่อไป (เรียกมันว่า x1) ที่น่าจะทำให้ f(x) ใกล้ศูนย์มากขึ้น ควรจะคำนวณอย่างนี้ครับ:

x1 = x0 – f(x0)/f'(x0) โดยที่ f'(x) คือ derivative ของ f(x) ครับ

ถ้าค่า x1 ทำให้ f(x) ไม่ใกล้ 0 พอ เราก็หา x2, x3, x4, … ไปเรื่อยๆจนเราพอใจว่าค่า f(xn) ใกล้ 0 พอแล้ว โดยที่ xn หาได้จาก xn-1 ดังนี้ครับ:

xn = xn-1 – f(xn-1)/f'(xn-1)

ถ้าจะทำการคำนวณด้วยวิธีของนิวตันใน Mathematica สามารถทำอย่างนี้ครับ:

newtonSolve[f_, guess_, steps_] := NestList[ #1 – f[#1]/f'[#1] &, guess, steps]

f คือฟังค์ชั่นที่เราจะหา f(x) = 0
guess คือค่าที่เราเดาตอนแรกว่า f(guess) น่าจะไม่ห่างจาก 0 นัก
steps คือจำนวนครั้งที่เราจะทำการทำวิธีของนิวตันซำ้ๆกัน

ยกตัวอย่างเช่น เราจะหาค่ารูทที่สองของสอง เราก็เขียนสมการ f(x) = x^2 -2 = 0 ก่อน เพราะค่า x เท่ากับรูทที่สองของสองจะแก้สมการนั้นพอดี:

f[x_] := x^2 – 2

แล้วเราก็เรียก newtonSolve ด้วยฟังค์ชั่น f โดยเดาค่า guess = 1 และให้ทำซ้ำสักห้าครั้ง:

newtonSolve[f, 1, 5]

แล้วเราก็ได้ผลดังนี้: {1, 3/2, 17/12, 577/408, 665857/470832, 886731088897/627013566048}

Mathematica ทำการคำนวณให้เป็นค่าเศษส่วนไม่มีทศนิยม เพราะเราเดาด้วยค่า 1 ซึ่งเป็นจำนวนที่ไม่มีการประมาณเข้ามาเกี่ยวข้อง ถ้าเราต้องการคำตอบเป็นเลขทศนิยม เราก็ควรเดาด้วยค่า 1.0 ดังนี้:

newtonSolve[f, 1.0, 5]

แล้วเราก็จะได้ผลดังนี้: {1., 1.5, 1.41667, 1.41422, 1.41421, 1.41421} ซึ่งเราจะเห็นว่า 1.41421 นั้นเป็นค่าประมาณของรูทที่สองของสองได้ดีทีเดียว

นักศึกษาลองไปทดลองดูครับ

มีใครอยากลองอธิบายว่า newtonSolve[f_, guess_, steps_] := NestList[ #1 – f[#1]/f'[#1] &, guess, steps] ทำงานอย่างไรใน Mathematica ไหมครับ เป็นการฝึกความเข้าใจเรื่อง Pure function และ Functional programming ครับ

******
education, math, Mathematica, ภาษาไทย

1, 11, 21, 1211, 111221 แล้วอะไรต่อไปเอ่ย

Leave a comment

เฉลยที่นี่

(เฉลยสั้นๆว่าตัวต่อไปให้อ่านตัวปัจจุบันดังๆว่ามีเลขอะไรกี่ตัวครับ เช่นตัวปัจจุบันคือ 1 ตัวต่อไปก็อ่านว่า หนึ่งหนึ่ง (มีหนึ่งหนี่งตัว) เขียนว่า 11
ตัวต่อไปอ่านว่า สองหนึ่ง (มีหนึ่งสองตัว) เขียนว่า 21
ตัวต่อไปอ่านว่า หนึ่งสองหนึ่งหนึ่ง (มีเลขสองหนึ่งตัว ตามด้วยเลขหนึ่งหนึ่งตัว) เขียนว่า 1211ตัวต่อไปอ่านว่า หนึ่งหนึ่ง หนึ่งสอง สองหนึ่ง (มีเลขหนึ่งหนึ่งตัว ตามด้วยเลขสองหนึ่งตัว ตามด้วยเลขหนึ่งสองตัว))

ผมสงสัยว่าจำนวนหลัก (หน่วย สิบ ร้อย พัน …) ในพจน์ที่ N จะใหญ่ประมาณ 1.324N นั่นก็คือพจน์ที่ N จะยาวขึ้นเร็วมากๆๆๆ