Computing an infinite series function symbolically in matlab

270 Views Asked by At

Let $u(x)$ be a function of $x$ defined as an infinite series (below). I would like to compute $u$ for a certain range of $x$.

Here are the values of the parameters, $\alpha=0.5$, $K_{\alpha} = 1$, $t = 10$, and $-40<x<40$.

$$u(x) = \frac{1}{\sqrt{4K_{\alpha}t^{\alpha}}} \sum_{k=0}^{\infty} \frac{(-1)^k}{k! \Gamma(1-0.5\alpha(k+1))} \left(\frac{x^2}{K_{\alpha}t^{\alpha}}\right)^{(k/2)}$$

Here is the matlab code which is used to compute $u$.

clear all
clc
format long e
syms x real % x represents spatial variable
syms xx real  
syms u_n real % u_n is going to contain the nth partial sum

n = 100; % number of terms used for the partial sum

T = 10; L = 40;
dx = 0.1; %x = (-L:dx:L);
x1 = (0:dx:L);
x2 = sort(-x1(2:length(x1)));
x3 = [x2 x1]; %J = length(x);


gma = 0.5; K_gma = 1; 

a1 = 1/sqrt(4*K_gma*T^gma);

gama = (-1).^((1:n)-1)./(factorial((1:n)-1).*gamma(1-0.5*gma*(1:n)));

xx = ( x.^2 / (K_gma*T^gma)).^(0.5*((1:n)-1));

u_n = a1*sum(gama.*xx);

U = subs(u_n, x3);

figure(1); 
plot(x3, U)   
axis([-10 10 0 0.35])
xlabel('x'), ylabel('u'),
title('\alpha = 0.5' ) 

My Problem is the following.

If you look at the entries of U they are not simplified, for example

U(1) = (13727875503963355228348566315747*(160*10^(1/2))^(5/2))/20769187434139310514121985316880384 - (12869883284965642519161564104547*(160*10^(1/2))^(1/2))/81129638414606681695789005144064 - (13945778289740548006536715712349*(160*10^(1/2))^(9/2))/42535295865117307932921825928971026432 + (2080489068865722820853778749805*(160*10^(1/2))^(13/2))/43556142965880123323311949751266331066368 - (8354591398111767502175328215655*(160*10^(1/2))^(17/2))/2854495385411919762116571938898990272765493248 + (17153085898078092321098620177869*(160*10^(1/2))^(21/2))/187072209578355573530071658587684226515959365500928 - (10182469541815633930394116379361*(160*10^(1/2))^(25/2))/6129982163463555433433388108601236734474956488734408704 + (7609449962283677024317653385533*(160*10^(1/2))^(29/2))/401734511064747568885490523085290650630550748445698208825344 - (7616888329793240347301396286771*(160*10^(1/2))^(33/2))/52656145834278593348959013841835216159447547700274555627155488768 + (10707430149749676615840652030947*(160*10^(1/2))^(37/2))/13803492693581127574869511724554050904902217944340773110325048447598592 - (10971265514290098609677073631071*(160*10^(1/2))^(41/2))/3618502788666131106986593281521497120414687020801267626233049500247285301248 + (16890130531936006833720592334403*(160*10^(1/2))^(45/2))/1897137590064188545819787018382342682267975428761855001222473056385648716020711424 - (5006656148415839444948741981661*(160*10^(1/2))^(49/2))/248661618204893321077691124073410420050228075398673858720231988446579748506266687766528 + (9337665196578709627158592793151*(160*10^(1/2))^(53/2))/260740604970814219042361048116400404614587954389239840081425977517360806369707098391474864128 - (13942885083686161163103989942067*(160*10^(1/2))^(57/2))/273406340597876490546562778389702670669146178861651554553221325801244124899921990402939147127881728 + (16926200182355256988797294071379*(160*10^(1/2))^(61/2))/286687326998758938951352611912760867599570623646035140467198604923365359511060601008752319138765710819328 - (16930333564573182737543504404317*(160*10^(1/2))^(65/2))/300613450595050653169853516389035139504087366260264943450533244356122755214669880763353471793250393988087676928 + (14117982217911386087503819557969*(160*10^(1/2))^(69/2))/315216049571155833698232320801148910440637914163723573343586347233965774171977684891314130039079325126453023922454528 - (19834846456000105249084985495733*(160*10^(1/2))^(73/2))/661055968790248598951915308032771039828404682964281219284648795274405791236311345825189210439715284847591212025023358304256 + (5923429014993948098540282457789*(160*10^(1/2))^(77/2))/346583711765101857447301773017885462929554634421977071896309947576827663475703202879996800763017447262173901370175446478621769728 - (6066539208105086020653891341079*(160*10^(1/2))^(81/2))/726838724295606890549323807888004534353641360687318060281490199180639288113397923326191050713763565560762521606266177933534601628614656 + (10734074814684009956198097625845*(160*10^(1/2))^(85/2))/3048582568667961163458591044719888970457615373696260889510895468384152088691177363398736428772941378085768487423248655171335913749304966119424 - (16518575770611103911004476082461*(160*10^(1/2))^(89/2))/12786682062094304179739022253232809188346257992355721833919106906625522642205759980012773798148063113870651109873281527379754908382364816614564560896 + (22246430892586454057441037263163*(160*10^(1/2))^(93/2))/53631231719770388398296099992823384509917463282369573510894245774887056120294187907207497192667613710760127432745944203415015531247786279785734596024336384 - (3296125192115571724961294089419*(160*10^(1/2))^(97/2))/28118211215894977392565865673037386617935606989386978956879722328823984879196799189494004288149317857187005691459505594520051662846839373056303219880407274094592 + (13912754125562980015419445309270293633256178869371469422876626746385997243379061407485*10^(1/2))/474284397516047136454946754595585670566993857190463750305618264096412179005177856 + 16717258617315192158890516822547512640982950817091919248532413798082315647077061526199/118571099379011784113736688648896417641748464297615937576404566024103044751294464

How can it be simplified. Thanks

1

There are 1 best solutions below

0
On BEST ANSWER

There are several ways to get a numerical value of a symbolic expression in Matlab.

  1. As commented by Ander Biguri, you can use double(U). For example, double(U(1)) returns 0.7706.
  2. The second solution is that you can use eval(U). For example, eval(U(1)) returns 0.7706.
  3. The third solution is that you can use vpa(U). For example, vpa(U(1)) returns 0.77057045019475968092311323544069.

    vpa can also be used with a specified precision. For example, vpa(U(1), 4) returns 0.7158
    vpa(U(1), 8) returns 0.77059168 vpa(U(1), 12) returns 0.770570447887

So, you can add the following line U = vpa(U,12); in your Matlab code after the line U = subs(u_n, x3);

It will solve your problem.