2
$\begingroup$

The background is that I'm solving a problem in numerical analysis and I seem to be halfway having solved the non-linear system and now getting to the part where I should apply solving differential equation.

https://math.stackexchange.com/questions/467164/how-to-proceed-solving-this-problem

The code so far is

z=[0:500:4000 5000:1000:12000];
data=[5050 4980 4930 4890 4870 4865 4860 4860 4865 4875 4885 4905 4920 4935 4950 4970 4990];
fun=@(p)(4800 + p(1))*ones(size(z)) +p(2)/1000*z+p(3)*exp(p(4)/1000*z)-data;
x0=[0 0 0 -1];
opt = optimset('MaxFunEvals',1000);
p=lsqnonlin(fun,x0,[],[],opt);
fitf=@(t)(4800 + p(1))*ones(size(t)) + p(2)/1000*t+ p(3)*exp(p(4)/1000*t);
tt=linspace(0,12000,1000);
plot(z,data,'r-',tt,fitf(tt),'b-');

And the system of equations is stated in the math question. Now I suppose that I should define a function something like

function f=fsys(x, z)
    f= ...

But first I suppose that I must rewrite the system from second order to first order which is just math theory, right? So I first rewrite the system to first order and then I define the function in matlab and finally something like

xend=...; z0=...
[X,Z]=ode45(@fsys,[0, xend],z0); zs=Z(end), plot(X,Z,'--')

Can you help me define fsys, tell me how to proceed, confirm or reject that I'm on the right way? Any hints if you don't want to help me but hint at how to learn matlab well enough to solve this "pretty simple math"?

I think I can transform the system in the usual way:

u = z'
u' = z''

The system becomes

u'=-q0*c'(z)/((c(z)^3)

Since z is a function of x, u and u' are also functions of x and thus this system till become a first-ordet system but how do I program the solution in matlab? How do I use that u=tan(b0)?

$\endgroup$
2
$\begingroup$

I take $p4=1$ and i get the other coefficients as $p_{1}=5.66$, $p_{2}=14.58$, $p_{3}=257.4$. So $c(z)$ and $c'(z)$ will be as

$c(z)=4800+5.66+(14.58)\frac{z}{1000}+(257.4)exp(-\frac{z}{1000})$ and

$c'(z)=(1287exp(\frac{z}{1000}))/5000+729/50000$

Let $z'=u$ and $z''=u'$ than we may write following differential eq.

$$ \frac{d^{2}z}{dx^{2}}=\left(\frac{c(2000)}{cos(7.8)}\right)^{2}\frac{c'(z)}{c(z)^{3}}, \qquad z(0)=2000,\qquad z'(0)=tan(7.8) $$

as

\begin{align*} z' & =u,\qquad z(0)=2000\\ u' & =\left(\frac{c(2000)}{cos(7.8)}\right)^{2}\frac{c'(z)}{c(z)^{3}},\qquad u(0)=tan(7.8) \end{align*}

So the codes for the problem should be as

    function dZ=sys(x,Z)
    c=@(z)4800 + 5.66 + (14.58)*z/1000+ (257.4)*exp(-z/1000); % c(z)
    c=c(2000);

% Z(1):=z
% Z(2):=u
    dZ=zeros(2,1);    % a column vector
    dZ(1)=Z(2);
    dZ(2)=(c/cosd(7.8))^2*((1287*exp(Z(1)/1000))/5000 + 729/50000)/...
        (4800 + 5.66 + (14.58)*Z(1)/1000+ (257.4)*exp(-Z(1)/1000))^3;
    end

We can call the above function by the following codes and plot the solution

>> x=0:0.5:3000;
>> [X,Z]=ode45(@sys,x,[2000 tand(7.8)]);
>> plot(X,Z(:,1),'r',X,Z(:,2),'b')  %Z(:,1) is z(x) and Z(:,2) is z'(x).

You see just z below enter image description here

$\endgroup$
  • $\begingroup$ Thanks! But we've gotten different coefficients. I have to check that the coefficients are correct. The matlab I run is z=[0:500:4000 5000:1000:12000]; data=[5050 4980 4930 4890 4870 4865 4860 4860 4865 4875 4885 4905 4920 4935 4950 4970 4990]; fun=@(p)(4800 + p(1))*ones(size(z)) +p(2)/1000*z+p(3)*exp(p(4)/1000*z)-data; x0=[0 0 0 -1]; opt = optimset('MaxFunEvals',1000); p=lsqnonlin(fun,x0,[],[],opt); fitf=@(t)(4800 + p(1))*ones(size(t)) + p(2)/1000*t+ p(3)*exp(p(4)/1000*t); tt=linspace(0,12000,1000); plot(z,data,'r-',tt,fitf(tt),'b-'); and then p(4) is not 1, is it something else? $\endgroup$ – Niklas Aug 31 '13 at 12:06
  • 1
    $\begingroup$ The coefficients i got would be more precise actually i didn't focus on coeffs compeletely, because you were trouble with conversion of the original ode to the system of ode so i focused on the ode45 usage. Anyway you can use your coefficients also. Just replace your coeffs in the codes. $\endgroup$ – Ömer Aug 31 '13 at 21:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.