> # The Existence and Uniqueness Theorem and Picard's Method # The existence and uniqueness theorem as a solution method # # (Recall: Given a function f(x,y) that is continuous in a rectangle # a # y'=f(x,y) # is the standard form, so here, > > > > # f(x,y)=x+y^2. > > # Picard's Method: # Here we limit ourselves to solutions through the point (0,0). That is, # we take x0=y0=0. If the starting point is different, a substitution # can fix that. Or, the method can be suitably modified. # Initial approximate solution: phi0(x):=f(0,0), a constant. # For this problem, the constant is zero=0+0^2. > # Rule to generate the next approximate solution: # phi_(x) is replaced with > x / | 2 | t + phi(t) dt | / 0 # > > > # Time Out...we need to do this by hand for a stretch: > > > > > # Time In: > # Functions in Maple: > g1:=x^2; 2 g1 := x > g2:=x->x^2; 2 g2 := x -> x > g3:=proc(x)x^2 end; g3 := proc(x) x^2 end > g4:=unapply(g1,x); 2 g4 := (x -> x) # TYPO Maple bug. It does the right thing, but prints it wrong. > g4(3); 9 > g4(5); 25 > g3(5); 25 > g2(5); 25 > g1(5); 2 x(5) > subs(x=5,g1); 25 > plot(g1,x=0..2); > plot(g2,0..2); # Now we need to get on to the problem at hand. How does Picard's method # work on our problem, y'=x+y^2? First, we need to tell Maple what # f(x,y) is. > > f:=(x,y)->x+y^2; 2 f := (x, y) -> x + y # Next, we need to tell Maple what the first iteration phi0 is. It's a # function, even if it's a constant function. > phi0:=x->0; phi0 := 0 # Next, we tell it the rule for Picard's method. This has been # specialized to the case y0=0. # # "Picard" is a Maple procedure (function) that takes as input a # FUNCTION, and returns a new, modified FUNCTION. The point is, that # each successive application of "Picard" gives a function which more # nearly satisfies the original differential equation. > picard:=proc(phi)unapply(int(f(t,phi(t)),t=0..x),x)end; picard := proc(phi) unapply(int(f(t, phi(t)), t = 0 .. x), x) end > > picard(phi0); 2 x -> 1/2 x > phi1:="; 2 phi1 := x -> 1/2 x > phi2:=picard(phi1); 2 5 phi2 := x -> 1/2 x + 1/20 x > phi3:=picard(phi2); 2 11 8 5 phi3 := x -> 1/2 x + 1/4400 x + 1/160 x + 1/20 x > phi4:=picard(phi3); 2 5 23 20 phi4 := x -> 1/2 x + 1/20 x + 1/445280000 x + 1/7040000 x 87 17 14 11 8 + -------- x + 3/49280 x + 7/8800 x + 1/160 x 23936000 > plot([phi0,phi1,phi2,phi3,phi4],0..1); > plot([phi0,phi1,phi2,phi3,phi4],1..1.2); # So, it appears that we have a sequence of functions which is # converging to `THE ANSWER'. # # What was all the fuss and bother in the blackboard presentation of the # theorem, about "there EXISTS" an interval in which the solution holds? # Why can't we go on like this indefinitely? > plot([phi0,phi1,phi2,phi3,phi4],1.2..1.5); > plot([phi0,phi1,phi2,phi3,phi4],1.5..1.8); > phi5:=picard(phi4):phi6:=picard(phi5); 2 5 11 8 phi6 := x -> 1/2 x + 1/20 x + 7/8800 x + 1/160 x 145237 23 1027 20 2169 17 14 + ------------ x + --------- x + --------- x + 1/9856 x 847813120000 670208000 167552000 30155989 29 68664301557182789 47 + ----------------- x + -------------------------------- x 16271845990400000 95793703809363897745408000000000 1608443 26 449004488539017 44 + -------------- x + ----------------------------- x 88172564480000 46724706256354436710400000000 3687821403309 41 2963358526269 38 + -------------------------- x + ------------------------- x 30340718348282101760000000 2035048181896970240000000 3106476321017 35 5142027483 32 + ------------------------ x + -------------------- x 187438648332615680000000 28638448943104000000 92 + 1/29563142485230611333120000000000000000 x 3043428839544096211423 71 + ----------------------------------------------- x 64608504706380752631259658211819520000000000000 142775831073384583 50 + ---------------------------------- x 2841961545239440444620800000000000 2550960910951100573 77 + ----------------------------------------------- x 54092865726474448609945994238361600000000000000 95 + 1/8249964449784667475148800000000000000000 x 57784903116340450781 65 + ------------------------------------------- x 2038235234045064590393112474419200000000000 1438295993 86 + ------------------------------------------- x 3726258133612793826126069760000000000000000 30663101735449147473 62 + ----------------------------------------- x 51533334171211497413356080332800000000000 454248291881180555815961 59 + -------------------------------------------- x 39820306023780651904501017686835200000000000 1689767509231341 83 + ----------------------------------------------- x 70184851861090053462629527060480000000000000000 398665529549292274449 74 + ------------------------------------------------ x 248982472865399136663961944771788800000000000000 17155880993855001653583 53 + ---------------------------------------- x 5206226868616528111114712514560000000000 194087104854719 80 + --------------------------------------------- x 164995243869738591742884577280000000000000000 287637 89 + ----------------------------------------- x 63787492792567322185302016000000000000000 93739057249414069429 56 + --------------------------------------- x 465462373159329654056119435264000000000 347320565202822647477 68 + --------------------------------------------- x 284217024590349862231818234757120000000000000 > plot([phi0,phi1,phi2,phi3,phi4,phi5,phi6],1.2..1.5); > plot([phi0,phi1,phi2,phi3,phi4,phi5,phi6],1.5..1.8); # Ouch. Even with all those terms, there is no clear limiting curve # emerging. > plot([phi0,phi1,phi2,phi3,phi4,phi5,phi6],1.8..2); # The fact of the matter is, the solution can shoot off to infinity # [That is, phi(0)=0, phi(x1)=Infinity] in finite time. The Picard # iterations, individually, will not. But they won't converge to a # limiting curve when x>=x1. Instead, the picture will look as it does # above. Each new iteration just produces a new curve, not particularly # close to the one that came before, and none of them any longer # signifying anything.