Allan Wittkopf

240 Reputation

7 Badges

18 years, 246 days

MaplePrimes Activity


These are replies submitted by Allan Wittkopf

Sorry, the plots got scrambled by html conversion, trying again:

  +              AAAAAAA                               AAAAAA                 
0.35            AA      AA                           AA      AA               
  +           AA          AA                       AA          AA             
0.3          A             AA                     AA            AA            
  +         A                A                   A                A           
  +        A                  A                 AA                AA          
0.25      A                   AA               A                    A         
  +      AA                    AA             AA                    AA        
  +      A                      A             A                      A        
0.2     A                        A           A                        A       
  +    A                          A         AA                         A      
  +    A                          A         A                          A      
0.15  A                            A       A                            A     
  +  AA                            AA     AA                            AA    
0.1  A                              A     A                              A    
  + AA                              AA   AA                              AA   
  + A                                A   A                                A   
0.05A                                A   A                                A   
  +A                                  A A                                  A  
  *                                   AAA                                   A 
  *--+--+--+--+--+--+---+--+--+--+--+--*--+--+--+--+--+---+--+--+--+--+--+--* 
0               0.5              1             1.5              2             

   +                                                                  AAAAAA  
   +                                                              AAAAA       
   +                                                           AAA            
0.5+                                                         AA               
   +                                                      AAA                 
   +                                                    AAA                   
0.4+                                                  AA                      
   +                                               AAA                        
   +                                            AAA                           
   +                                        AAAA                              
0.3+                               AAAAAAAAA                                  
   +                           AAAA                                           
   +                        AAA                                               
0.2+                     AAA                                                  
   +                   AA                                                     
   +                AAA                                                       
   +              AAA                                                         
0.1+            AA                                                            
   +         AAA                                                              
   +    AAAAA                                                                 
  -******+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+- 
                 0.5             1             1.5             2              
Hello Ian, Yes, this is a bug, and it has to do with the fact that the dependency lists of the two functions are different. Consider the following equivalent problem: > dsys1 := [diff(u(y,x),y)=0,diff(u(y,x),x,x)=0,v(y,x)=0]: > DEtools[initialdata](dsys1); table([Finite = [u(y[0], x[0]) = _C1, D[2](u)(y[0], x[0]) = _C2], Infinite = [] ]) > dsys2 := [diff(u(x,y),y)=0,diff(u(x,y),x,x)=0,v(x,y)=0]: > DEtools[initialdata](dsys2); table([Finite = [u(x[0], y[0]) = _C1, D[1](u)(x[0], y[0]) = _C2], Infinite = [] ]) So as you see, if the dependency lists of 'u' and 'v' are the same, the problem does not appear. A fix is being worked on.
I am unsure what you are referring to in this post. The desired outcome, is that a return value? A displayed value? How are 'b','A','C' assigned? Are you assigning 'a'? Could you please provide a code snippet that illustrates what you are doing/trying to do, and we may be able to provide more assistance.
By any chance are either of the 'x' or 't' in the x(t) being used in the dsolve call local variables? Consider this: > f := proc() local x; x; end proc: > val := f(); val := x > val - x; x - x Note that the 'x' declared as local is a different 'x' than just 'x' at the top level. If this is the case, you either need to return the 'x' along with the dsolve solution in your procedure, and used that in the call to odeplot, or declare the 'x' (and the 't') to be global variables, and take some care that they are not assigned.
What difference schemes to use entirely depend on how your domain is defined, and what type of boundary conditions you have. This is especially true in a problem as you have shown, where you have mixed order coupled equations. For example, if your domain is semi-infinite, you will probably want to use one-sided difference equations (and probably a transform), while for a finite region, you will probably want to use central differences. Even resticting to finite problems, it is still not so clear as to what you need. If for the third equation, you know S11 and S11_x2 for the left-most boundary, then you will want to use a one-sided difference, but if you know S11 on both the right and left boundaries, you instead want to use a central difference. I'm not even discussing the issue of whether the discretization is stable, but in some cases, even though central differences seem to be best, they can lead to bad solutions, so this is also a concern. So what is your domain, and what boundary conditions do you have? - Allan
Hello Thomas, Sorry, I should have been more clear. You need to pass all the B[i,j] used by your 'W' routine explicitly as individual parameters. So if N was 3, you would need to pass: i,k,j,l,B[1,1],B[1,2],B[1,3],B[2,1],.... or the 4 integer parameters, and the 9 B[i,j] parameters that form the system. When the numerical DE solver needs to evaluate the equations, it does so by passing in an array of strictly numerical values, and getting out an array of strictly numerical values - it cannot handle passing of arrays of the dependent variables to other procedures. Now you can explicitly code a general procedure to handle this type of problem by writing Maple code for the evaluation of the rhs of the DE system. In your case, here is what I would do (I will generally provide examples with N=3). I have attached the worksheet content below:

1) Determine an ordering of the B[i,j] that maps these to a single sequence

of variables - I would just use the same sequence order that you have

provided with your double sequence call, essentially row-major order.

> N := 3:

> dvars := [seq(seq(B[i,k](t),k=1..N),i=1..N)];

Maple Equation

2) Now code a dsolve/numeric evaluation procedure explicitly to do

   what you want:

> dproc := proc(n,t,Y,YP)

> local N,W,B,i,j,k,l;

>

>    # System is n=N^2, so determine N from n

>    for N from 1 while N^2

>    if N^2<>n then

>       error "unexpected number of variables";

>    end if;

>

>    # Now transfer values to the B[i,j] for convenience

>    B := hfarray(1..N,1..N):

>    for i to N do

>       for j to N do

>          B[i,j] := Y[N*(i-1)+j];

>       end do;

>    end do;

>

>    # Now declare and compute 'W'

>    W := hfarray(1..N,1..N,1..N,1..N);

>    # Some complicated computation for W[i,j,k,l] here that can depend on

>    # the B[i,j], on 't' on N, etc.

>    # For now, use your example, W=1/N^2

>    for i to N do

>       for j to N do

>          for k to N do

>             for l to N do

>                W[i,j,k,l] := 1/N^2;

>             end do;

>          end do;

>       end do;

>    end do;

>

>    # Now compute the value for the derivative of B[i,j] from the B[i,j],

>    # W[i,j,k,l], N, and t placing the values in the output array YP

>    for i to N do

>       for k to N do

>          YP[N*(i-1)+k] := add(add(W[j,l,i,k]*B[j,l]

>                                  -W[i,k,j,l]*B[i,k],l=1..N),j=1..N);

>       end do;

>    end do;

> end proc:

3) Now at this point, you have a mapping of the variables, and a procedure

that can be used to compute the diff(B[i,j](t),t) for the i,j with the

W you program, so all we need now is initial conditions, and to perform

the call to generate the dsolve/numeric procedure.

> ini := Array([seq(seq(i-j,j=1..N),i=1..N)]);

> dsn := dsolve(numeric,procedure=dproc,initial=ini,start=0,procvars=dvars);

Maple Equation

Maple Equation

And the procedure can be used to compute the solutions:

> dsn(0);

Maple Equation

> dsn(1);

Maple Equation

> Maple Equation

This post generated using the online HTML conversion tool
Download the original worksheet

Hello Thomas, Sorry, I should have been more clear. You need to pass all the B[i,j] used by your 'W' routine explicitly as individual parameters. So if N was 3, you would need to pass: i,k,j,l,B[1,1],B[1,2],B[1,3],B[2,1],.... or the 4 integer parameters, and the 9 B[i,j] parameters that form the system. When the numerical DE solver needs to evaluate the equations, it does so by passing in an array of strictly numerical values, and getting out an array of strictly numerical values - it cannot handle passing of arrays of the dependent variables to other procedures. Now you can explicitly code a general procedure to handle this type of problem by writing Maple code for the evaluation of the rhs of the DE system. In your case, here is what I would do (I will generally provide examples with N=3). I have attached the worksheet content below:

1) Determine an ordering of the B[i,j] that maps these to a single sequence

of variables - I would just use the same sequence order that you have

provided with your double sequence call, essentially row-major order.

> N := 3:

> dvars := [seq(seq(B[i,k](t),k=1..N),i=1..N)];

Maple Equation

2) Now code a dsolve/numeric evaluation procedure explicitly to do

   what you want:

> dproc := proc(n,t,Y,YP)

> local N,W,B,i,j,k,l;

>

>    # System is n=N^2, so determine N from n

>    for N from 1 while N^2

>    if N^2<>n then

>       error "unexpected number of variables";

>    end if;

>

>    # Now transfer values to the B[i,j] for convenience

>    B := hfarray(1..N,1..N):

>    for i to N do

>       for j to N do

>          B[i,j] := Y[N*(i-1)+j];

>       end do;

>    end do;

>

>    # Now declare and compute 'W'

>    W := hfarray(1..N,1..N,1..N,1..N);

>    # Some complicated computation for W[i,j,k,l] here that can depend on

>    # the B[i,j], on 't' on N, etc.

>    # For now, use your example, W=1/N^2

>    for i to N do

>       for j to N do

>          for k to N do

>             for l to N do

>                W[i,j,k,l] := 1/N^2;

>             end do;

>          end do;

>       end do;

>    end do;

>

>    # Now compute the value for the derivative of B[i,j] from the B[i,j],

>    # W[i,j,k,l], N, and t placing the values in the output array YP

>    for i to N do

>       for k to N do

>          YP[N*(i-1)+k] := add(add(W[j,l,i,k]*B[j,l]

>                                  -W[i,k,j,l]*B[i,k],l=1..N),j=1..N);

>       end do;

>    end do;

> end proc:

3) Now at this point, you have a mapping of the variables, and a procedure

that can be used to compute the diff(B[i,j](t),t) for the i,j with the

W you program, so all we need now is initial conditions, and to perform

the call to generate the dsolve/numeric procedure.

> ini := Array([seq(seq(i-j,j=1..N),i=1..N)]);

> dsn := dsolve(numeric,procedure=dproc,initial=ini,start=0,procvars=dvars);

Maple Equation

Maple Equation

And the procedure can be used to compute the solutions:

> dsn(0);

Maple Equation

> dsn(1);

Maple Equation

> Maple Equation

This post generated using the online HTML conversion tool
Download the original worksheet

It may then be more suitable to use a piecewise approximation to the solution - this allows it to be used more flexbly later on. So for your problem, generate the piecewise solution for a given set of parameters, and z=0..2. As a simple example: > sol := op([2,2],dsolve({diff(y(x),x)=-y(x),y(0)=1},numeric, > range=0..2,output=piecewise)); sol := { { undefined , x < 0. 0.999621690653975 - 0.972619698891573 x 2 + 0.486310338843125 (x - 0.0277621231571284) 3 - 0.162115462190059 (x - 0.0277621231571284) 4 + 0.0402113832410770 (x - 0.0277621231571284) , ... Is this suitable?
It may then be more suitable to use a piecewise approximation to the solution - this allows it to be used more flexbly later on. So for your problem, generate the piecewise solution for a given set of parameters, and z=0..2. As a simple example: > sol := op([2,2],dsolve({diff(y(x),x)=-y(x),y(0)=1},numeric, > range=0..2,output=piecewise)); sol := { { undefined , x < 0. 0.999621690653975 - 0.972619698891573 x 2 + 0.486310338843125 (x - 0.0277621231571284) 3 - 0.162115462190059 (x - 0.0277621231571284) 4 + 0.0402113832410770 (x - 0.0277621231571284) , ... Is this suitable?
To answer each of the questions in turn: (1) why cannot I evaluate the value of f(t) at t=5 by using f(5)? You just need to do this a little differently. Consider: > A := 100; v := 50; rho := 0.996; r := 0.9; theta := 0.993; A := 100 v := 50 rho := 0.996 r := 0.9 theta := 0.993 > x :=5000; c := 200; m := 100000; l := 5000; x := 5000 c := 200 m := 100000 l := 5000 > solution := dsolve([ > diff(s(t),t) = A - A * rho^((1 - r^(theta * t)) * x) - v, > diff(f(t),t)=((c - s(t)) / l) * m + x, > diff(h(t),t) = x, > f(0) = 0, s(0) = 0, h(0) = 0], > numeric, output=operator); bytes used=8307108, alloc=917336, time=0.41 solution := [t = (proc(t) ... end proc), f = (proc(t) ... end proc), h = (proc(t) ... end proc), s = (proc(t) ... end proc)] > f := eval(f,solution); f := proc(t) ... end proc So 'f' is now a procedure that can be used to evaluate the 'f' values at different points: > f(100); 7 -0.399982377673172718 10 > f(1000); 9 -0.489993379968400420 10 Please be warned though that it is not such a good idea to assign 'f' then use it again. Generally it is better to use something else like 'fsol' for 'f'. Note that if you try to apply dsolve again your original system of DEs is no longer a system of DEs, but is a system with a 'f' function, and attempting to specify the initial condition for 'f' gives an identity: > f(0) = 0; 0. = 0 (2) For each fixed x, there are curves s(t), f(t), and h(t). Given s(tau) = c, I want to find f(tau)=?. How can I do that? Do I need to find tau first, then find f(tau)? How? Yes, continuing on from prior, we have: > s := eval(s,solution); s := proc(t) ... end proc > fsolve(s(tau)=100); -0.6237217782 > val := fsolve(s(tau)=100,tau=0..100); val := 3.000783529 > f(val); 24992.3520950843448 So use 'fsolve' to compute the tau value at which s(tau)=100, then evaluate with 'f'. (3) My ultimate goal is for x in a range, say x=500..5000, find the minimum of f(tau) (since for each x there is a f(tau)). How can I do that? The easiest way to do this is to delay the evaluation of the DE solution, leaving 'x' a variable, then create an interface procedure to evaluate it at a specific time value, or as follows: > A := 100; v := 50; rho := 0.996; r := 0.9; theta := 0.993; A := 100 v := 50 rho := 0.996 r := 0.9 theta := 0.993 > c := 200; m := 100000; l := 5000; c := 200 m := 100000 l := 5000 > solution := dsolve([ > diff(s(t),t) = A - A * rho^((1 - r^(theta * t)) * x) - v, > diff(f(t),t)=((c - s(t)) / l) * m + x, > diff(h(t),t) = x, > f(0) = 0, s(0) = 0, h(0) = 0], > numeric, output=operator); solution := [t = (proc(t) ... end proc), f = (proc(t) ... end proc), h = (proc(t) ... end proc), s = (proc(t) ... end proc)] > f := eval(f,solution); f := proc(t) ... end proc Now the 'f' solution cannot be evaluated until 'x' is assigned: > f(1); Error, (in f) global 'x' must be assigned to a numeric value before obtaining a solution Now the interface procedure, say you want values at t=100, then do: > tval := 100: > getval := proc(xval) global x; > if type(xval,'numeric') then > x := xval; f(tval); > else > 'getval'(xval); > end if; > end proc: Now 'getval' is a function that computes f(t=100;x), and it can be plotted: > plot(getval,500..5000); * -2.4e+06 *A -2.6e+06 +A -2.8e+06 + A + A -3e+06 + A -3.2e+06 + AA -3.4e+06 + AA -3.6e+06AA + AA -3.8e+06 AA + AA + AAA -4e+06 AAAAAA AAAAAAAAAAAAAAAAAAAAA +-+--+--+--+---+--+-**********************************+---+--+--+--+---+--+ 1000 2000 3000 4000 5000 Or solved: > x := 'x': > Optimization[Minimize](getval(x),x=500..5000); 7 [-0.413681014780686982 10 , [2480.90973324573588 = 2480.90973324573588]] - Allan
1 2 3 4 5 6 Page 6 of 6