แบบจำลองการแพร่กระจายเชื้อมาลาเรียในร่างกายผู้ป่วย

แบบจำลองคณิตศาสตร์สำหรับการแพร่กระจายเชื้อมาลาเรียในร่างกายผู้ป่วยนั้นมีผู้เสนออยู่หลายแบบครับแต่แบบที่ผมคิดว่ามันเจ๋งที่สุดก็คือแบบจำลองที่คิดโดย White NJ, Chapman D, Watt G. The effects of multiplication and synchronicity on the vascular distribution of parasites in falciparum malaria. Trans R Soc Trop Med Hyg. 1992;86(6):590-7. PubMed PMID: 1287908. ครับ

ตัวอย่างข้างล่างนี้ผมลองเอาแบบจำลองมาเปรียบเทียบกับข้อมูลจริงที่ได้จากผู้ป่วยครับ

ข้อมูลจากผู้ป่วยนี้ก็เป็นอะไรที่น่าสนใจครับ เพราะได้มาจากการเฝ้านับเชื้อมาลาเรียจากผู้ป่วยที่เป็นซิฟิลิสครับ ในสมัยก่อนผู้ป่วยที่เป็นซิฟิลิสจะถูกรักษาด้วยการทำให้ติดเชื้อมาลาเรียครับ คนที่คิดวิธีนี้ก็ได้รางวัลโนเบลด้วยนะครับ (https://en.wikipedia.org/wiki/Julius_Wagner-Jauregg) สมัยนี้คงทำไม่ได้แล้ว เสี่ยงตายมากๆ

มาใช้ Mathematica ในการเรียนการสอนกัน ด้วยคำสั่ง Manipulate 3

จากตอนที่แล้ว เดี๋ยวผมจะพูดถึงการใช้งาน slider กับตัวควบคุมอื่นๆ ของ Manipulate ครับ

ดูตัวอย่างไปก่อนนะครับ

Manipulate[
Module[{anglegraph, maingraph},
anglegraph[th_, showtext_] := Show[
Graphics[{
{Lighter[Gray, 0.5], Circle[{0, 0}, 1]},
{Darker[Green, 0.2], Thick, Circle[{0, 0}, 1, {0, th}]},
{Lighter[Gray, 0.5], Line[{{0, 0}, {Cos[th], Sin[th]}}]},
{Red, Thick, Line[{{Cos[th], 0}, {Cos[th], Sin[th]}}]},
If[showtext,
Rotate[Text[
Style[N[Sin[th]], 11], {Cos[th] – 0.15, Sin[th]/2}],
3 Pi/2, {Cos[th] – 0.15, Sin[th]/2}], {}],
{Lighter[Red, 0.3], Dashing[Medium],
Line[{{Cos[th], Sin[th]}, {2, Sin[th]}}]}
}],
PlotRange -> 1, ImageSize -> 145, BaseStyle -> {12}, Axes -> True,
Ticks -> {{-1, 1}, {-1, 1}}, PlotRange -> {{-1, 1}, {-1, 1}},
PlotRangePadding -> 0.25];

maingraph[th_, showtext_] := Module[{},
Show[Plot[{Sin[x]}, {x, 0.0001, th},
PlotRange -> {{0, 2 Pi}, {-1, 1}},
PlotRangePadding -> {0, 0.25},
ImagePadding -> {{30, 12}, {0, 0}}, PlotRangeClipping -> False,
PlotStyle -> Gray,
Ticks -> {Table[{n Pi/4, n Pi/4}, {n, 0, 8}],
Table[n, {n, -1, 1, 1/2}]},
GridLines -> {Table[{n Pi/4, Lighter[Gray, 0.7]}, {n, -2, 8}],
Table[{n, Lighter[Gray, 0.7]}, {n, -1, 1, 1/2}]},
ImageSize -> {Automatic, 145}],
Graphics[{If[showtext,
Rotate[Text[Style[N[Sin[th]], 11], {th + 0.15, Sin[th]/2}],
3 Pi/2, {th + 0.15, Sin[th]/2}], {}],
{Lighter[Red, 0.3], Dashing[Medium],
Line[{{th, Sin[th]}, {-2, Sin[th]}}]},
{Darker[Green, 0.2], Thick, Line[{{0, 0}, {th, 0}}]},
{Red, Thick, Line[{{th, 0}, {th, Sin[th]}}]}
}], AspectRatio -> Automatic, BaseStyle -> {12}]];

DynamicModule[{pt = {Cos[ptctrl], Sin[ptctrl]}, pt2 = {ptctrl, 0}},
Labeled[Grid[{
{LocatorPane[Dynamic[pt,
{(pt = {Cos[pt2[[1]]], Sin[pt2[[1]]]}) &,
(pt = Normalize[#];
pt2 = {If[pt2 == {2 Pi, 0}, 2 Pi,
Mod[ArcTan[#[[1]], #[[2]]], 2 Pi]], 0}) &,
(pt = Normalize[#]; ptctrl = pt2[[1]]) &}],
Dynamic[
anglegraph[
If[pt2 == {2 Pi, 0}, 2 Pi,
Mod[ArcTan[pt[[1]], pt[[2]]], 2 Pi]], showvalue]]],
LocatorPane[Dynamic[pt2,
{(pt2 = {If[pt2 == {2 Pi, 0}, 2 Pi,
Mod[ArcTan[pt[[1]], pt[[2]]], 2 Pi]], 0}) &,
(pt2 = {#[[1]], 0}; pt = {Cos[#[[1]]], Sin[#[[1]]]}) &,
(pt2 = {#[[1]], 0}; ptctrl = #[[1]]) &}],
Dynamic[
maingraph[
If[pt2 == {2 Pi, 0}, 2 Pi,
Mod[ArcTan[pt[[1]], pt[[2]]], 2 Pi]], showvalue]]]}},
Spacings ->
0], {Row[{Style[“Illustrating “, “Label”, 20, Gray],
Text@Style[“sin(“, Red, 24],
Text@Style[“x”, Italic, Bold, Darker[Green, 0.3], 24],
Style[“)”, Red, 24]}],
Style[“”, 10, Lighter[Gray, 0.7], “Label”]}, {{Top,
Left}, {Bottom, Right}}]]
],
{{showvalue, False, “show value”}, {False, True}},
{{ptctrl, Pi/6, “angle”}, 0, 2 Pi},
TrackedSymbols :> {showvalue, ptctrl}]

ส่วนอันนี้เป็นตัวอย่างที่เราสามารถลากเขียนกราฟได้ครับ

Manipulate[
points = If[Last[points] == p, points, points~Join~{p}];
ListLinePlot[points, PlotRange -> {{-100, 100}, {-100, 100}},
PlotLabel -> (“จำนวนจุด ” <> ToString@Length[points])],
Button[“Clear”, clear = True; points = {{0, 0}};
p = {0, 0};], {{p, {0, 0}}, Locator},
Initialization :> {
points = {{0, 0}};
clear = False;
}
]

 

Manipulate[Quiet@DynamicModule[{g=9.81,m=1.,deqns,xsol,ysol,times,periods,periodNumber},

deqns={(*parabolic*){a (1+x[t]^2) Derivative[2][x][t]==-\[Mu] Derivative[1][x][t]/m-x[t] (g+2 a Derivative[1][x][t]^2)},(*elliptic*){\[Mu] Derivative[1][x][t]/m+Derivative[2][x][t]==x[t] (-g Sqrt[1-x[t]^2]/a+Derivative[1][x][t]^2/(-1+x[t]^2))}};
{xsol,times}=Reap@NDSolveValue[{deqns[[i]],x[0]==.999,Derivative[1][x][0]==0,WhenEvent[x[t]==0,Sow[t]]},x,{t,tMax}];
ysol[t_]:={(*parabolic path*)a (xsol[t]^2-1),(*elliptic path*)-a Sqrt[1-xsol[t]^2]}[[i]];
periodNumber=Position[times[[1]],Nearest[times[[1]],tt][[1]]][[1,1]]-1;
periods=Differences[times[[1]]];
Column[{Row[{“swing# “,periodNumber,”: “,If[periodNumber!=0,periods[[periodNumber]],0],” seconds”}],Dynamic@Show[{Plot[{a (x^2-1),-a Sqrt[1-x^2]}[[i]],{x,-1,1},PlotStyle->Thick],Graphics[{Red,Disk[.98 {xsol[tt],ysol[tt]},.02]}]},AspectRatio->Automatic,ImageSize->400,PlotRange->{{-1,1},{-1.25,0.05}}]}]],{{i,1,”path”},{1->”parabolic”,2->”elliptic”},Setter},{{a,1},.5,1.25,.001,Appearance->”Labeled”},{{\[Mu],.15},0,.25,.01,Appearance->”Labeled”},Control[{{tt,0},0,.5 tMax,.01,Trigger,AnimationRate->1}],{{tMax,100},None}]

ตัวอย่างการใช้ Mathematica แก้ปัญหาแบบง่ายๆ

พอดีมีคนถามว่า
จะหาจำนวนนับที่เล็กที่สุด ที่มีคุณสมบัติครบทั้ง 2 ข้อต่อไปนี้ได้ยังไงครับ???
I
(1) ขึ้นต้นด้วยเลข 1
(2) เมื่อสลับตัวเลขหลักแรก (ซึ่งก็คือ 1) กับตัวเลขหลักสุดท้าย แล้วจำนวนใหม่มีค่าเป็น 3 เท่าของจำนวนเดิม
II
(1) ขึ้นต้นด้วยเลข 1
(2) เมื่อย้ายตัวเลขหลักแรก (ซึ่งก็คือ 1) ไปต่อหลังตัวเลขหลักสุดท้าย แล้วจำนวนใหม่มีค่าเป็น 3 เท่าของจำนวนเดิม

(http://www.pantip.com/cafe/wahkor/topic/X11362457/X11362457.html)

อันนี้เป็นคำสั่ง Mathematica แบบง่ายๆ ที่ลองเขียนเพื่อแก้ปัญหานี้
I
swapInt1[tmp_Integer] := FromDigits@ ReplacePart[IntegerDigits[tmp], {1 -> Last@IntegerDigits[tmp], Length@IntegerDigits[tmp] -> First@IntegerDigits[tmp]}]
Catch[If[swapInt1@# == 3*#, Throw[#]] & /@ Range[1, 1000000, 1]]
คำตอบที่ได้จะเป็น Null ทั้งหมด นั่นหมายความว่ามันไม่มีคำตอบ
II
swapInt2[tmp__Integer] := FromDigits@RotateLeft@IntegerDigits@tmp

Catch[If[swapInt2@# == 3*#, Throw[#]] & /@ Range[1, 1000000, 1]]

คำตอบที่ได้คือ 142857
————————————————

จาก http://www.pantip.com/cafe/wahkor/topic/X11435196/X11435196.html

เราเรียก 6 28 496 8128….
ว่าเป็นจำนวนสมบูรณ์
เพราะมันมีคุณสมบัติที่น่าสนใจคือ
6     มีตัวประกอบแท้คือ 1,2,3
6  = 1+2+3
28   มีตัวประกอบแท้คือ 1,2,4,7,14
28  = 1+2+4+7+14 28 = 1+2+3+4+5+6+7
496 มีตัวประกอบแท้คือ 1,2,4,8,16,31,62,124,248
496 = 1+2+4+8+16+31+62+124+248
496 = 1+2+3+…+31

โดยตัวเลขพวกนี้ มันมีนิยามว่า มันอยู่ในรูป 2n-1(2n-1) เมื่อ n เป็นจำนวนนับที่ทำให้  2n-1 เป็นจำนวนเฉพาะ ซึ่งจากนิยามนี้เองเราสามารถเขียนคำสั่ง Mathematica ง่ายๆเพื่อหาตัวเลขพวกนี้ได้ว่า

(# (# + 1))/2 & /@ Select[2^Range[100] – 1, PrimeQ]

(#(#+1))/2 นี้ได้มาจากการให้ (2n-1) = x –> 2n =x+1  ดังนั้น  2n-1= (x+1)/2

——————————————
จาก
http://www.pantip.com/cafe/wahkor/topic/X11494743/X11494743.html
http://www.pantip.com/cafe/wahkor/topic/X11495190/X11495190.html

มีเลขอยู่หกหลักที่ไม่เหมือนกันซักตัว
เลขชุดนี้สามารถถอดรากที่สองได้เป็นจำนวนเต็ม
เมื่อกลับเลขชุดนี้ ก็จะยังถอดรากที่สองได้เป็นจำนวนเต็มเหมือนเดิม
เช่น MINTLA —> ถอดรากได้เป็นจำนวนเต็ม       ALTNIM —> ถอดรากได้เป็นจำนวนเต็ม
ถามว่า มีจำนวนใดบ้างที่เข้ากฏนี้และมีกี่จำนวน
ตัวอย่าง Mathematica code สำหรับปัญหานี้
(* ตั้งแต่ 3 – 10 หลัก แบบไม่ซ้ำ *)
For[i = 3, i <= 10, i++,
Print[“==”, i “==”];
num = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
test = Developer`ToPackedArray@ Select[Permutations[num, {i}], #[[1]] != 0 &];
For[j = 1, j <= Length@test, j++,
If[FractionalPart@N@Sqrt@FromDigits@test[[j]] == 0 &&
FractionalPart@N@Sqrt@FromDigits@Reverse@test[[j]] == 0,
Print[FromDigits@test[[j]]]]
] ]

==3 ==
169
961
==4 ==
1089
9801
==5 ==
12769
96721
==6 ==
==7 ==
1238769
9678321
==8 ==
==9 ==
==10==

—————————–
http://www.pantip.com/cafe/wahkor/topic/X11514071/X11514071.html

ABCDEF แทนจำนวนใดได้บ้าง?
กำหนดให้
(1) ABCDEF แทนจำนวนซึ่งมี 6 หลัก
(2) ตัวอักษรแต่ละตัว แทนเลขโดด 1 – 9 ซึ่งอาจซ้ำกันได้
(3) AB ถอดรากที่สองแล้วได้จำนวนเต็ม
(4) CD ถอดรากที่สองแล้วได้จำนวนเต็ม
(5) EF ถอดรากที่สองแล้วได้จำนวนเต็ม
(6) ABCDEF ถอดรากที่สองแล้วได้จำนวนเต็ม

ตัวอย่างคำสั่ง Mathematica สำหรับปัญหานี้ครับ
num = {1, 2, 3, 4, 5, 6, 7, 8, 9};
test = Tuples[num, 6];

For[j = 1, j <= Length@test, j++,
If[FractionalPart@N@Sqrt@FromDigits@test[[j, # ;; # + 1]] & /@ {1, 3, 5} == {0, 0, 0} &&
FractionalPart@N@Sqrt@FromDigits@test[[j]] == 0,
Print[FromDigits@test[[j]]]]
]

166464
646416

——–

จาก http://www.pantip.com/cafe/wahkor/topic/X12638167/X12638167.html

 

ทำด้วย Mathematica

NestList[Total[IntegerDigits[#]^3] &, ตัวเลขที่หารด้วย 3 ลงตัว, จำนวนรอบจากข้อ (2-4) ]

เช่น เริ่มที่เลข 3 จำนวน 10 รอบ NestList[Total[IntegerDigits[#]^3]&, 3, 10]

{3, 27, 351, 153, 153, 153, 153, 153, 153, 153, 153}

หรือ เริ่มที่เลข 123345 (ตามตัวอย่าง) 10 รอบ NestList[Total[IntegerDigits[#]^3]&, 123345, 10]

{123345, 252, 141, 66, 432, 99, 1458, 702, 351, 153, 153}

มันมาหยุดอยู่ที่ 153

 

P. falciparum in a patient during treatment with artesunate

หลังจากที่ดูของคนอื่นมานาน ก็เลยลองทำส่งบ้าง 🙂

มันเป็นโมเดลที่จำลองว่าปริมาณของเชื้อมาลาเรียในร่างกายของคนไข้จะเปลี่ยนไปอย่างไรระหว่างการรักษาด้วยยาอาทีซูเนท

http://demonstrations.wolfram.com/AModelOfPlasmodiumFalciparumPopulationDynamicsInAPatientDuri/

มาใช้ Mathematica ในการเรียนการสอนกัน ด้วยคำสั่ง Manipulate 2

ต่อจากครั้งที่แล้วครับ(http://www.sakngoi.com/?p=143)

คำสั่ง Manipulate นี้ Mathematica เพิ่มเข้ามาตั้งแต่ version 6 เป็นคำสั่งที่สามารถทำให้เราสามารถดูได้ครับว่ารูปแบบของกราฟหรือค่าจากการคำนวณอะไรบางอย่างที่เราสนใจจะเปลี่ยนไปอย่างไรถ้าหากเราเพิ่มหรือลดค่าที่เราสนใจที่เป็นส่วนหนึ่งในการคำนวณนั้น

เช่น อยากรู้ว่ากราฟของ sin\left(\omega\theta\right) ที่plot ตั้งแต่ \theta เท่ากับ 0-2\pi จะเปลี่ยนไปอย่างไรถ้าค่า \omega ค่อยเพิ่มขึ้นจาก 1-10 หรือ อย่างปัญหาในฟิสิกส์อยากรู้ว่าจะต้องยิงวัตถุจากจุดยอดของทรงกลมรัศมีขนาดหนึ่งด้วยมุมและอัตราเร็วเท่าใดวัตถุนั้นจึงจะเฉียดผิวของทรงกลมนี้พอดี และปัญหาอื่นๆอีกมากมาย สามารถดูตัวอย่างได้ที่ http://demonstrations.wolfram.com/

รูปแบบของคำสั่ง Manipulate ก็ตามที่แสดงด้านล่างนี้ครับ

รูปแบบคำสั่งนี้เป็นแบบ “จัดเต็ม”ครับ 🙂  แต่เราก็สามารถใส่แบบสั้นๆได้ครับ เช่นในการกำหนดตัวแปร อยากจะใส่เพียง {ตัวแปร,ค่าน้อยสุด,ค่ามากสุด} ก็ได้ครับ step ก็เป็นค่าอย่างเช่นจาก 1 ถึง 10 เราจะให้ค่าค่อยเพิ่มขึ้นครั้งล่ะเท่าไหร่จาก 1 ไปจนถึง 10  ส่วน options ก็จะเกี่ยวข้องกับการแสดงผลที่เกี่ยวข้องกับ slider หรือตัวควบคุมอื่นๆ เช่น Animator, Checkbox, ColorSetter, ColorSlider, InputField, Manipulator, PopupMenu, RadioButton หรือ RadioButtonBar, Setter หรือ SetterBar, Slider2D, Trigger and VerticalSlider ซึ่งเดี๋ยวจะพูดถึงทีหลังครับ

อันนี้เป็นตัวอย่างของกราฟ sin\left(\omega\theta\right)  โดย Manipulate จะสร้าง slider สำหรับการเปลี่ยนค่า \omega มาให้

Manipulate[

Plot[Sin[\omega \theta],{\theta,0,2\pi}]

,{\omega,1,10}

]

 

ส่วนอันนี้ก็เป็นการใช้ PolarPlot กับ cos(9 \theta) โดย \theta ค่อยๆเพิ่มจาก 0.01 ถึง \pi

Manipulate[
PolarPlot[Cos[9 \theta], {\theta, 0, T},
PlotRange -> {{-1, 1}, {-1, 1}}],
{{T, 0.01, “\theta (radian)”}, 0.01, \pi, Appearance -> “Labeled”}]

หรือจะลอง plot กราฟของฟังชั่นทางตรีโกนมิติ เช่น

Manipulate[
Plot[amp fun[freq x], {x, 0, 10}, PlotRange -> {-3, 3},
PlotStyle -> color, PlotLabel -> fun], {freq, 1, 5}, {amp, 1,
5}, {fun, {Sin, Cos, Tan, Csc, Sec, Cot}}, {{color,
Red}, {Purple -> “Purple”, Green -> “Green”, Blue -> “Blue”, Yellow -> “Yellow”}}]

ส่วนอันนี้ก็เป็นการเอาไปประยุกต์กับปัญหาฟิสิกส์ที่ว่าจะต้องยิงวัตถุจากจุดยอดของทรงกลมรัศมีขนาดหนึ่งด้วยมุมและอัตราเร็วเท่าใดวัตถุนั้นจึงจะเฉียดผิวของทรงกลมนี้พอดี(http://mpec.sc.mahidol.ac.th/forums/index.php/topic,345.0.html)

Manipulate[
Show[{Graphics[Circle[{0, 0}, r, {0, Pi}], Axes -> True],
Plot[Tan[theta] x – (9.8/(2 u^2 Cos[theta]^2)) x^2 + r, {x, 0,
r + 10}, PlotStyle -> Red,
PlotRange -> {{-r – 10, r + 10}, {0, r + 10}}]}],
{{u, 4.95, “initial speed(m/s)”}, 0.01, 50, 0.0001, Appearance -> “Labeled”},
{{theta, 0.5236, “launch angle(radian)”}, 0, Pi/2, 0.0001, Appearance -> “Labeled”},
{{r, 5, “circle radius(m)”}, 1, 10, Appearance -> “Labeled”}]

ส่วนอันนี้เป็น demonstration project ที่ผมลองทำส่งไปที่เวบ Wolfram ครับ

http://demonstrations.wolfram.com/AModelOfPlasmodiumFalciparumPopulationDynamicsInAPatientDuri/

ต่อ..สร้างโปรแกรมด้วย ManipulateMaker