Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab Code

  1. Sep 14, 2015 #1
    Hi, I have to solve one equation on Matlab. The equation is: (p-c)*psi'-psi=0 I have to find p. Where psi is:


    and psi' is the derivative of psi. I've written this code:
    Code (Matlab M):

    rho = 0.04; %In the image is lambda
    c = 25;
    eta = 0.02;
    mu = 40;
    sigma  = 0.57;
    alpha = -rho/eta;

    syms p t
    F = t^(-alpha-1)*exp(-t^2/2+((p-mu)/sigma*sqrt(2*eta))*t); %Funzione integranda non derivata
    F_d = t^(-alpha-1)*(t*sqrt(2*eta))/sigma*exp(-t^2/2+((p-mu)/sigma*sqrt(2*eta))*t); %Funzione integranda derivata
    D = (exp(((p-mu)/sigma*sqrt(2*eta))^2/4)/gamma(-alpha))*int(F,t,0,inf); %Funzione cilindrica non derivata
    D_d = (exp((((p-mu)/sigma*sqrt(2*eta))^2)/4)/gamma(-alpha))*int(F_d,t,0,inf)+((eta*(p-mu))/(4*(sigma)^2))*(exp((((p-mu)/sigma*sqrt(2*eta))^2)/4)/gamma(-alpha))*int(F,t,0,inf); %Funzione cilindrica derivata
    psi = exp((eta*(p-mu)^2)/(2*(sigma)^2))*D;
    psi_d = (eta*(p-mu))/(sigma^2)*exp((eta*(p-mu)^2)/(2*(sigma)^2))*D + exp((eta*(p-mu)^2)/(2*(sigma)^2))*D_d; %Derivata di psi
    Theta = (p-c)*psi_d-psi;

    FreeBoundary = solve(Theta,p)
    I've calculated the derivatives by myself, and I hope they are ok. What do you think about the code? It works only if rho>eta. Moreover for different values of sigma it give me some strange results.
    Do you think the code is right and the derivatives are correct?
  2. jcsd
  3. Sep 19, 2015 #2
    Thanks for the post! This is an automated courtesy bump. Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post?
  4. Oct 22, 2015 #3


    User Avatar
    Gold Member

    I think the problem with this post is the question posed:

    "Do you think the code is right and the derivatives are correct?"

    It is unreasonable to request that someone manually check all of your work and then debug your code on this scale. This could easily take a few hours, which is why nobody has responded.

    Let's separate the math from the code: if you have a problem with the derivative or equations, post that in the Calculus forums. Once you nail down the expression you want to solve, I suggest you repost here with some updated code. The code you posted didn't run for me, so there is at least a little more work to be done:

    Code (Text):

    Warning: Cannot solve symbolically. Returning a numeric
    approximation instead.
    > In solve at 306

    FreeBoundary =

    Any updates or other information you can provide?
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook