taylorode(f,x,y,x0,y0,n) := block([yy,xx,k:1,listarith:true], yy : taylor(y0,x,x0,1), xx : taylor(x,x,x0,1), while k < n do ( yy: y0 + integrate(subst(cons(x=xx,map("=",y,yy)),f),x,x0,x), yy:taylor(yy,x,x0,k), k : k + 1, xx:taylor(x,x,x0,k)), map("=",y,yy)); keepfloat : true; /* MY PROBLEM INPUT DATA (for my program) WAS diff ( y , x , 1 ) = expt ( 2.0 * x + 1.0 , sin ( x ) ) * ( cos ( x ) * ln ( 2.0* x + 1.0 )+ ( 2.0 * sin ( x ) ) / ( 2.0 * x + 1.0 ) ); ! Digits := 32; max_terms := 30; ! x_start := 0.1; x_end := 1.0 ; diff(y,0,exact_soln_y(x_start)); glob_look_poles := true; glob_max_iter := 1000000; glob_max_h := 0.00001; ! exact_soln_y := proc(x) return(expt(2.0*x+1.0,sin(x))); end; */ /* This is my problem using taylor and integrate with 30 terms*/ elapsed_real_time(); h(x) := ((2.0*x+1.0) ^ (sin(x))); elapsed_real_time(); taylorode([ ( 2.0 * x + 1.0) ^( sin ( x ) ) * ( cos ( x ) * log ( 2.0* x + 1.0 )+ ( 2.0 * sin ( x ) ) / ( 2.0 * x + 1.0 )) ],x,[y],0.1,[(2.0*0.1+1.0)^(sin(0.1))],30); elapsed_real_time(); /* I didn't know how to turn the returned value into a function so I edited the output and added it to the program. Hint welcome */ g(x) := block([], 1.01836844605944 + .3541314832678168 * (x - 0.1) + 1.601687062117677 * (x - 0.1) ^2 - .8547087765566383 * (x - 0.1) ^ 3 + 2.271776575775109 * (x - 0.1) ^ 4 - 5.134287865791374 * (x - 0.1) ^ 5 + 13.97222222222217 * (x - 0.1) ^ 6 - 26.01666666666458 * (x - 0.1) ^ 7 + 48.57420634914174 * (x - 0.1) ^ 8 - 91.45515872848526 * (x - 0.1) ^ 9 + 173.1665663206988 * (x - 0.1) ^ 10 - 329.4505338889471 * (x - 0.1) ^ 11 + 629.361372314316 * (x - 0.1) ^ 12- 1206.495472393729 * (x - 0.1) ^ 13 + 2319.86503960143 * (x - 0.1) ^ 14 - 4472.409704978651 * (x - 0.1) ^ 15 + 8642.03292878619 * (x - 0.1) ^ 16 - 16731.62708448044 * (x - 0.1) ^ 17 + 32440.9032281923 * (x - 0.1) ^ 18 - 62924.60292145748 * (x - 0.1) ^ 19 + 121791.6712878852 * (x - 0.1) ^ 20 - 233892.2314154419 * (x - 0.1) ^ 21 + 440759.2422171746 * (x - 0.1) ^ 22 - 800045.7726936718 * (x - 0.1) ^ 23 + 1361553.762834966 * (x - 0.1) ^ 24 - 2097241.582954726 * (x - 0.1) ^ 25 + 2798845.538938754 * (x - 0.1) ^ 26 - 3062586.070419776 * (x - 0.1) ^ 27 + 2545421.348878588 * (x - 0.1) ^ 28 - 1415068.435589704 * (x - 0.1) ^ 29 ); elapsed_real_time(); printf(true,"g is the taylorode series ~%"); printf(true,"h is the closed form solution ~%"); elapsed_real_time(); g(0.1); elapsed_real_time(); h(0.1); elapsed_real_time(); /* The following are like taking my step size */ g(0.10001); elapsed_real_time(); h(0.10001); elapsed_real_time(); /* This is impressive but I need to cause the output of taylorode to be a function, and also a way to get ahold of the individual coefficients. I am sure there are simple ways. I tried a couple ways but failed. The output was: Maxima 5.31.3 http://maxima.sourceforge.net using Lisp CLISP 2.49 (2010-07-07) Distributed under the GNU Public License. See the file COPYING. Dedicated to the memory of William Schelter. The function bug_report() provides bug reporting information. (%i1) batch("to2.max") read and interpret file: /home/dennis/mastersource/mine/taylorode/to2.max (%i2) taylorode(f,x,y,x0,y0,n):=block([yy,xx,k:1,listarith:true], yy:taylor(y0,x,x0,1),xx:taylor(x,x,x0,1), while k < n do (yy:integrate(subst(cons(x = xx,map("=",y,yy)),f),x,x0, x) +y0,yy:taylor(yy,x,x0,k),k:1+k,xx:taylor(x,x,x0,k)), map("=",y,yy)) (%o2) taylorode(f, x, y, x0, y0, n) := block([yy, xx, k : 1, listarith : true], yy : taylor(y0, x, x0, 1), xx : taylor(x, x, x0, 1), while k < n do (yy : integrate(subst(cons(x = xx, map("=", y, yy)), f), x, x0, x) + y0, yy : taylor(yy, x, x0, k), k : 1 + k, xx : taylor(x, x, x0, k)), map("=", y, yy)) (%i3) keepfloat:true (%o3) true (%i4) elapsed_real_time() (%o4) 0.362975 (%i5) h(x):=(1.0+2.0*x)^sin(x) sin(x) (%o5) h(x) := (1.0 + 2.0 x) (%i6) elapsed_real_time() (%o6) 0.495445 (%i7) taylorode([(1.0+2.0*x)^sin(x)*(2.0*sin(x)/(1.0+2.0*x) +cos(x)*log(1.0+2.0*x))],x,[y],0.1, [(1.0+2.0*0.1)^sin(0.1)],30) (%o7)/T/ [y = 1.01836844605944 + .3541314832678168 (x - 0.1) 2 3 + 1.601687062117677 (x - 0.1) - .8547087765566383 (x - 0.1) 4 5 + 2.271776575775109 (x - 0.1) - 5.134287865791374 (x - 0.1) 6 7 + 13.97222222222217 (x - 0.1) - 26.01666666666458 (x - 0.1) 8 9 + 48.57420634914174 (x - 0.1) - 91.45515872848526 (x - 0.1) 10 11 + 173.1665663206988 (x - 0.1) - 329.4505338889471 (x - 0.1) 12 13 + 629.361372314316 (x - 0.1) - 1206.495472393729 (x - 0.1) 14 15 + 2319.86503960143 (x - 0.1) - 4472.409704978651 (x - 0.1) 16 17 + 8642.03292878619 (x - 0.1) - 16731.62708448044 (x - 0.1) 18 19 + 32440.9032281923 (x - 0.1) - 62924.60292145748 (x - 0.1) 20 21 + 121791.6712878852 (x - 0.1) - 233892.2314154419 (x - 0.1) 22 23 + 440759.2422171746 (x - 0.1) - 800045.7726936718 (x - 0.1) 24 25 + 1361553.762834966 (x - 0.1) - 2097241.582954726 (x - 0.1) 26 27 + 2798845.538938754 (x - 0.1) - 3062586.070419776 (x - 0.1) 28 29 + 2545421.348878588 (x - 0.1) - 1415068.435589704 (x - 0.1) + . . .] (%i8) elapsed_real_time() (%o8) 5.59718 (%i9) g(x):=block([], -1415068.435589704*(x-0.1)^29+2545421.348878588*(x-0.1)^28 -3062586.070419776*(x-0.1)^27 +2798845.538938754*(x-0.1)^26 -2097241.582954726*(x-0.1)^25 +1361553.762834966*(x-0.1)^24 -800045.7726936718*(x-0.1)^23 +440759.2422171746*(x-0.1)^22 -233892.2314154419*(x-0.1)^21 +121791.6712878852*(x-0.1)^20 -62924.60292145748*(x-0.1)^19 +32440.9032281923*(x-0.1)^18 -16731.62708448044*(x-0.1)^17 +8642.03292878619*(x-0.1)^16 -4472.40970497865*(x-0.1)^15 +2319.86503960143*(x-0.1)^14 -1206.495472393729*(x-0.1)^13 +629.361372314316*(x-0.1)^12 -329.4505338889471*(x-0.1)^11 +173.1665663206988*(x-0.1)^10 -91.45515872848526*(x-0.1)^9 +48.57420634914174*(x-0.1)^8 -26.01666666666458*(x-0.1)^7 +13.97222222222217*(x-0.1)^6 -5.134287865791374*(x-0.1)^5 +2.271776575775109*(x-0.1)^4 -.8547087765566383*(x-0.1)^3 +1.601687062117677*(x-0.1)^2 +.3541314832678168*(x-0.1) +1.01836844605944) 29 (%o9) g(x) := block([], - 1415068.435589704 (x - 0.1) 28 27 + 2545421.348878588 (x - 0.1) - 3062586.070419776 (x - 0.1) 26 25 + 2798845.538938754 (x - 0.1) - 2097241.582954726 (x - 0.1) 24 23 + 1361553.762834966 (x - 0.1) - 800045.7726936718 (x - 0.1) 22 21 + 440759.2422171746 (x - 0.1) - 233892.2314154419 (x - 0.1) 20 19 + 121791.6712878852 (x - 0.1) - 62924.60292145748 (x - 0.1) 18 17 + 32440.9032281923 (x - 0.1) - 16731.62708448044 (x - 0.1) 16 15 + 8642.03292878619 (x - 0.1) - 4472.40970497865 (x - 0.1) 14 13 + 2319.86503960143 (x - 0.1) - 1206.495472393729 (x - 0.1) 12 11 + 629.361372314316 (x - 0.1) - 329.4505338889471 (x - 0.1) 10 9 + 173.1665663206988 (x - 0.1) - 91.45515872848526 (x - 0.1) 8 7 + 48.57420634914174 (x - 0.1) - 26.01666666666458 (x - 0.1) 6 5 + 13.97222222222217 (x - 0.1) - 5.134287865791374 (x - 0.1) 4 3 + 2.271776575775109 (x - 0.1) - .8547087765566383 (x - 0.1) 2 + 1.601687062117677 (x - 0.1) + .3541314832678168 (x - 0.1) + 1.01836844605944) (%i10) elapsed_real_time() (%o10) 5.689776 (%i11) printf(true,"g is the taylorode series ~%") g is the taylorode series (%o11) false (%i12) printf(true,"h is the closed form solution ~%") h is the closed form solution (%o12) false (%i13) elapsed_real_time() (%o13) 6.261538 (%i14) g(0.1) (%o14) 1.01836844605944 (%i15) elapsed_real_time() (%o15) 6.345724 (%i16) h(0.1) (%o16) 1.01836844605944 (%i17) elapsed_real_time() (%o17) 6.4296 (%i18) g(0.10001) (%o18) 1.018371987534441 (%i19) elapsed_real_time() (%o19) 6.577728 (%i20) h(0.10001) (%o20) 1.0183719881022 (%i21) elapsed_real_time() (%o21) 6.70559 (%o22) to2.max */