#C12.txt, March 3, 2014. Picard's proof of the existence #of solutions to 1st order diffeq Help:=proc(): print(`DS(f,x,y,x0,y0)`): print(`Picard(f,x,y,x0,y0,N), Z1(f,x,y,cur) `): print(`Z(f,x,y,y0,N) , ERC(f,x) `): end: #The cheating way (using Maple) to solve the diff.eq #y'(x)=f(x,y(x)) y(x0)=y0 DS:=proc(f,x,y,x0,y0) local s: s:=dsolve({diff(y(x),x)=subs(y=y(x),f), y(x0)=y0},y(x)): if s=NULL then RETURN(FAIL): else RETURN(op(2,s)): fi: end: #Picard(f,x,y,x0,y0,N): the first N+1 terms #in the Picard sequence that is supposed to #converge to the solution of y'(x)=f(x,y(x)), y(x0)=y0 Picard:=proc(f,x,y,x0,y0,N) local L,n,yn,t: L:=[y0]: yn:=y0: for n from 1 to N do yn:=y0+int( subs({x=t,y=subs(x=t,yn) },f), t=x0..x): L:=[op(L),yn]: od: L: end: #Z1(f,x,y,cur): if cur is the current Maclaurin #polynomial approximating y'(x)=f(x,y(x)) finds #the poly appx. of one degree higher #in the Picard sequence that is supposed to #converge to the solution of y'(x)=f(x,y(x)), y(x0)=y0 #Right now f is a polynomial of x and y #Version of March 7, 2014 correcting semantical bugs spotted by #Andrew Lohr Z1:=proc(f,x,y,cur) local n,cur1,c: if cur=0 then n:=0: else n:=degree(cur,x): fi: cur1:=cur+c*x^(n+1): c:=solve(coeff(diff(cur1,x)-subs(y=cur1 ,f),x,n)=0,c): cur+c*x^(n+1): end: #Z(f,x,y,y0,N): the first N+1 MacLaurin polynomials #n x of sol. of y'(x)=f(x,y(x)), y(0)=y0 #Version of March 7, 2014 correcting semnatical bugs spotted by #Andrew Lohr Z:=proc(f,x,y,y0,N) local L,cur,n: if diff(f,y)=0 then print(`the function f only depends on x, so the solution is the indefinite integral`): RETURN(y0+int(f,x=0..x)): fi: L:=[y0]: cur:=y0: for n from 1 to N do cur:=Z1(f,x,y,cur): L:=[op(L),cur]: od: L: end: #given a MacLaurin polynomial f of x #appx. its radius of convergence emprically #with an error given ERC:=proc(f,x) local n: n:=degree(f,x): evalf([coeff(f,x,n-1)/coeff(f,x,n), coeff(f,x,n-2)/coeff(f,x,n-1)]): end: